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FOREWORD 


In September and October of 1985 Vice President Zhang Chi of the Chinese 
Aeronautical Establishment led a delegation of management and technical spe- 
cialists to the United States. Arrangements for the delegation were supervised 
by Lavonne Parker of NASA. Included In the activities of this delegation was 
the NASA - Chinese Aeronautical Establishment (CAE) Symposium on Combustion 
Fundamentals at the NASA Lewis Research Center In Cleveland, Ohio, from 
September 23 - 27, 1985. This symposium was one of several exchanges resulting 
from an agreement, signed on January 31, 1979, between the United States of 
America and the People's Republic of China to exchange technology In the area 
of civil aeronautics. The technical sessions of this symposium were conducted 
under the joint chairmanship of Mr. Zhou Xlaoqlng of the CAE and Dr. Edward J. 
Mularz of NASA. Following the presentation of overviews of combustion research 
activities In each of the organizations, 12 technical papers were given by com- 
bustion specialists from both CAE and NASA. Small group discussions were then 
held to clarify and further examine the presented material. These technical 
presentations are bound In this publication as a record of the proceedings of 
this joint symposium. 

In addition to the symposium participants, a number of technical observers 
were present. For the Chinese, these were the remaining members of the CAE 
delegation. Other observers were combustion specialists from Industry and 
local universities and some members of the Lewis technical staff. Following 
the symposium, the Chinese delegation toured combustion and related facilities 
at a number of relevant Industries and university establishments. 

Technical and other exchanges were open and fruitful with the technical 
experts using some similar methods, Instrumentation, and equipment. Common 
correlations and approaches and similar academic background existed for a 
number of the technical experts resulting In an apparent mutual understanding 
of the technology presented by both the CAE and the NASA combustion experts. 


Melvin J. Hartmann 
NASA Lewis Research Center 
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COMBUSTION RESEARCH IN THE INTERNAL FLUID MECHANICS DIVISION 


Edward J. Mularz 
Propulsion Directorate 

U.S. Army Aviation Research and Technology Actlvlty-AVSCOM 
NASA Lewis Research Center 
Cleveland, Ohio 


At the NASA Lewis Research Center, combustion research Is being conducted 
In the Internal Fluid Mechanics Division. The research Is organized Into three 
main functions focusing on the fluid dynamics related to aeropropulslon systems 
(fig. 1). The first function, computational methods, looks at Improved algo- 
rithms and new computational fluid dynamic techniques to solve Internal flow 
problems, Including heat transfer and chemical reactions. This area also looks 
at using expert systems and parallel processing as they might be applied to 
solving Internal flow problems. 

The second function Is fundamental experiments. These experiments can 
generate benchmark data In support of computational models and numerical codes, 
or they can focus on the physical phenomena of Interest to obtain a better 
understanding of the physics or chemistry Involved as a preamble to models and 
computer codes. 

Computational applications Is the third function In the Internal Fluid 
Mechanics Division. New flow codes are validated against available experi- 
mental data, and they are used as a tool to Investigate performance of real 
engine hardware. Since the geometry may be quite complex for the system being 
analyzed, large grids and computer storage may be required. The hardware of 
Interest Includes combustion chambers, high-speed Inlets, and turbomachinery 
components (e.g., a centrifugal compressor). 

The goal of this research Is to bring computational fluid dynamics to a 
state of practical application for the aircraft engine Industry. As shown In 
figure 2, the approach Is to have a strongly Integrated computational and 
experimental program for all the disciplines associated with the gas turbine 
and other aeropropulslon systems by advancing the understanding of flow 
physics, heat transfer, and combustion processes. The computational and 
experimental research Is Integrated In the following way: the experiments that 

are performed provide an empirical data set so that physical models can be 
formulated to describe the processes that are occurring - for example, turbu- 
lence or chemical reaction. These experiments also form a data base for those 
who are doing code development by providing experimental data against which 
the codes can be verified and assessed. Models are generated as closure to 
some of the numerical codes, and they also provide physical Insight for exper- 
iments. At the same time, codes which solve the complete Navler-Stokes 
equations can be used as a kind of numerical experiment from which far more 
extensive data can be obtained than ever could be obtained experimentally. 

This could provide physical Insight Into the complex processes that are taking 
place. These codes are also exercised against experimental data to assess the 
accuracy and applicability of models (e.g., the turbulence model). We feel 
that a fully Integrated computational-experimental research program Is more 
productive than other approaches and that It Is the most desirable way of 
pursuing our goal . 
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Figure 3 Is a cutaway view of a hypothetical combustor which Illustrates 
the typical complex fluid mechanics and combustion features. The flows are 
highly three-dimensional with turbulence levels, In many cases, comparable In 
magnitude to the bulk velocity. Liquid fuels are Injected as a spray which 
then undergoes vaporization and mixing. The chemical reaction which occurs 
causes changes In density and fluid mechanics properties and can cause the 
formation of a solid phase (soot) with Its attendant hlgh-radlatlon heat 
transfer. An understanding of these physical processes Is needed before 
accurate numerical codes can be built and used as a predictive tool In the 
design process. In addition, the numerical methods for three-dimensional 
flows need Improvements In accuracy and efficiency In order to properly 
simulate the features of these flows. 

This then Is the framework of the combustion research program. As shown 
In figure 4, the program Is divided Into four elements: advanced numerics, 
fuel sprays, fluid mixing, and radiation chemistry. The research In each of 
these elements Is focused on the long-range objective of developing numerical 
codes that can be used as a predictive tool to describe both two- and three- 
dimensional flows. Once the codes can be used with confidence. Industry will 
be able to Integrate them Into their combustor design system. 

Before discussing the specific research currently taking place, I would 
like to describe the major computing hardware now being used at NASA Lewis to 
compute reacting Internal flows. Lewis currently has a CRAY 1-S computer which 
Is tied to an IBM 370 system for Input/output. It has a relatively small 
storage capacity, but a very high calculation speed. In the near future a 
CRAY X-MP will be Installed (fig. 5). It will replace the CRAY 1-S and will 
provide a modest Increase In computer speed. Core storage, however, will 
Increase dramatically from 2 million words to 36 million words, and a high- 
speed solid-state device, which Is practically equivalent to core memory, will 
also be Included. In addition, the new equipment will have two processors 
which will give Lewis some ability to start using algorithms that can take 
advantage of simultaneous processing and thus Increase the computation speed 
even further. 


In addition to the computers located at Lewis, NASA Is building a "super 
computer" at Its Ames Research Center called the Numerical Aerodynamic 
Simulator (NAS) (fig. 6). The NAS computer will be accessed at both NASA 
Langley and NASA Lewis through satellite linkage with a UNIX operating system. 
In Its current configuration, the NAS computer Is a CRAY 2 with 256 megawords 
of memory, which will enable It to make very large calculations of chemically 
reacting Internal fluid flows. It will have four processors. In addition, as 
the state-of-the-art of computers advances changes will be made In the NAS 
computer to reflect these advances. 


The current research In combustion can be summarized In six major activ- 
ities. The first activity Is called Improved numerical methods for complex 
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This numerical error must be reduced before the accuracy of physical models can 
be assessed. The grids used In calculations that are performed with state-of- 
the-art codes are very sensitive to the flow direction relative to that grid, 
and at high angles this sensitivity leads to serious numerical errors. Tech- 
niques to reduce numerical diffusion, or error, and to virtually eliminate It 
from calculations are being pursued. Work Is also underway to establish 
second-order accurate closure models for turbulent reacting flow. 
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The second activity Involves the development of techniques for making 
predictive calculations of chemically reacting flow. One of the more promising 
techniques, which Is really In Its Infancy, Is called direct numerical simu- 
lations, or ONS. In this technique, the Navler-Stokes equations are solved 
directly without any modeling of the turbulence. This technique Is currently 
being applied to reacting shear layers. Although this Is a relatively simple 
flow, there Is nevertheless a lot of complexity associated with It. Much work 
remains to be done with this technique, but the results to date have been very 
promising. 

Another technique that has been used In several kinds of flows Is called 
the random vortex method. This technique accounts for the vortlclty generated 
at the wall of the confined flow and solves the vortlclty equation without any 
turbulence closure modeling. Calculated results of flow over a rearward 
facing step show many of the characteristics seen In high-speed movies of 
turbulent reacting flow experiments. 

The third activity Involves benchmark experiments for code development and 
verification. Two-phase flow research Is currently underway. Detailed data 
are required for code assessment, and Instruments which Lewis has helped to 
develop now show much promise of being able to make the appropriate measure- 
ments. (Those Instruments are also being used to support Icing research.) 
Research Is also being conducted on numerical calculation of two-phase flow. 
With the data from the experiments, an assessment of the current code capa- 
bility will be made to guide future code development efforts. In addition, an 
extensive set of experiments Is being conducted to look at the mixing of 
dilution jets Into a cross stream In a channel. A substantial range of para- 
metric variables has been studied, and a very complete set of data has been 
established. 

The next activity Is computer code applications. Currently available 
codes are applied to real systems. The study of flow In the transition sec- 
tion of a reverse flow combustor Is an example of such work. Here a flow In 
which fuel has already been burned has to undergo a 180° annular turn. Very 
strong secondary flows arise, and the analysis Is very complex. Another 
example Is the application of the Lewis-developed general chemical kinetics 
programs and chemical equilibrium codes to problems of practical Interest. 
These codes are used throughout the world. 

The next activity Is code development for thermochemical properties and 
kinetic rates. Thermodynamic calculations for real gases are being performed 
In support of future combustion models. The modeling of chemical kinetic com- 
putations for complex reactions Is also being actively developed In support of 
future computer codes. 

The final activity In the current combustion research program Is to 
advance the understanding of chemical mechanisms In reacting flow. Chemical 
kinetic rates are measured using a shock tube facility. The data from this 
experiment are used to assess the kinetic models for various fuels of Interest 
In aeropropulslon systems. Detailed characteristics of one- and two-dimen- 
sional controlled flames are also being established. We are Interested In 
studying the behavior of some of the minor species of flames and In looking at 
the soot nucleatlon growth and eventual soot consumption In flames. Radiation 
heat transfer Is dominant In combustors and other aeropropulslon systems. The 
control of soot nucleatlon and growth Is essential If radiation heat transfer 
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Is to be reduced. Also, some of the minor intermediates of combustion are 
being measured In an attempt to understand the detailed physics Involved In 
flame fronts. 
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APPLICATIONS 

Figure 1. - Fundamental aeropropulsion computational and experimental research in 
the internal fluid mechanics division. 



CLOSURE 


Figure 2. - Integrated computational-experimental methodology for propulsion 
aerothermodynamics research. Objective - advance the understanding of 
flow physics, heat transfer, and combustion processes which are funda- 
mental to aeropropulsion. 
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Figure 3. - Illustration of the typical flow phenomena in a gas turbine combustor. 
Typical flow is fully three-dimensional, has high turbulence levels, has chem- 
ical reactions and heat release, and occurs in two phases with vaporization. 
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Figure 4 - Fundamental combustion research plan. 
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Figure 5. - Comparison of Cray 1-S and Cray X-MP computer capabilities. 


NUMERICAL AERODYNAMIC SIMULATION 



Figure 6.- Numerical aerodynamic simulator (NAS) computer capabilities and Lewis remote workstation. 
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N87-20269 

NUMERICAL STUDY OF COMBUSTION PROCESSES IN AFTERBURNERS 

Zhou Xlaoqlng 

Chinese Aeronautical Establishment 
Beijing, People's Republic of China 

and 

Zhang Xlaochun 

Shenyang Aeroengine Research Institute 
Shenyang, People's Republic of China 


Mathematical models and numerical methods are presented for computer 
modeling of aeroengine afterburners. A computer code GEMCHIP Is described 
briefly. The algorithms SIMPLER, for gas flow predictions, and DROPLET, for 
droplet flow calculations, are Incorporated In this code. The block correction 
technique Is adopted to facilitate convergence. The method of handling 
Irregular shapes of combustors and flameholders Is described. The predicted 
results for a low-bypass-ratio turbofan afterburner In the cases of gaseous 
combustion and multiphase spray combustion are provided and analyzed, and 
engineering guides for afterburner optimization are presented. 


INTRODUCTION 

Thus far the design and development of air-breathing engine combustors 
have been based mainly on ad hoc tests. But the past 10 years have seen a 
boom of numerical fluid mechanics and combustion. The computer modeling tech- 
nique Is beginning to establish Itself In territory which has been dominated 
exclusively by the empirical technique. It Is believed that the combination 
of numerical modeling with the semi-empirical method and the modern diagnostic 
technology will greatly promote combustion study and combustor design and 
development. 

The objective of this paper (which Is part of the authors' continuing 
effort to make computer models of multiphase turbulent combustion processes In 
combustors and furnaces (refs. 1 and 2 and private communication with X. Zhang 
and H.H. Chiu, University of Illinois at Chicago, 1984)) Is to present mathe- 
matical models and numerical methods for combustor flows and to apply these 
methods to the computer prediction of aeroengine afterburner characteristics. 

The long-ignored droplet-turbulent-dlffuslon model plays an Important role 
In droplet dispersion and species distribution. It has now been Incorporated 
In this modeling. The K-e turbulence model Is modified and extended to account 
for multiphase turbulence effects. A hybrid, finite-rate chemical-reaction 
model based on a global Arrhenius law and a turbulent mixing-rate model Is used 
judiciously to predict the combustion rate of multiphase premixed turbulent 
flames. A system of conservation equations of the Eulerlan type Is used for 
both gaseous flow and multlslzed droplet flow predictions. To solve these 
equations, the authors developed a computer code called GEMCHIP (general, 
elliptic-type, multiphase, combustion-heat-transfer, and Interdiffusion 
program), In which the SIMPLER algorithm (ref. 3) for gaseous flow prediction 
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and the DROPLET (private communication with X. Zhou and H.H. Chiu, University 
of Illinois at Chicago, 1983) procedure for droplet flow calculation were 
Incorporated. The block correction technique and the alternating-direction 
reclprocatlng-sweep line-by-line TDMA (tri-diagonal matrix algorithm) method 
was adopted to facilitate convergence of the Iteration processes. 

The afterburner under study Is a low-bypass-ratio turbofan engine augmen- 
tor of the axlsymmetrlc type with a bluff body flameholder. The method of 
simulating the Irregularly shaped flow domain and the flameholder configuration 
Is described briefly and the Inlet and boundary conditions are presented. The 
usual wall functions (ref. 4) are used to bridge the near wall regions where 
local Reynolds numbers are very low. 

The numerical results of gas flow, droplet flow fields, and the spray 
flame structures are presented. A study of the parametric sensitivity of the 
combustion efficiency Is made for both gas combustion and multiphase combus- 
tion. Engineering guides are given for afterburner design optimization. 


SYMBOLS 

Cp specific heat of gas at constant pressure 

E activation energy 

h total enthalpy 

J work-heat equivalence 

K turbulent kinetic energy 

L latent heat of vaporization 

P gas pressure 

q heat value of fuel 

R u universal gas constant 

rj, droplet radius 

T temperature 

U,V velocities In X and R directions, respectively 
W molecular weight 

Y species mass concentration 

r gamma function; transport coefficient 

Y stoichiometric coefficient 
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I turbulence dissipation rate 

\ heat conductivity 

v viscosity 

p density 

§ mixture fraction 

Nondlmenslonal numbers: 

D] first Oamkohler number 

Ec Eckert number 

Eu Eulerlan number 

G c spray group combustion number 

G(j spray group aerodynamic drag number 

G v spray group droplet preheating number 

L m |< mass fraction of liquid In the mixture 

Nu Nusselt number 

Pe Pec let number 

Pr Prandtl number 

Re Reynolds number 

Sc Schmidt number 

Typ ratio of flow residence time to turbulence dissipation time 
c p ratio of turbulent kinetic energy to mean flow value 

a turbulent Schmldt/Prandtl number 

Subscripts: 
b boiling point 

eff effective 

F fuel 

g gas 

k Kth size group of droplets 

8, liquid 

9 



M,m mean 

o initial or Inlet condition 

ox oxygen 

p product 

T,t turbulent; temperature 


MATHEMATICAL FORMULATION 

The Eulerlan scheme Is used to construct elliptic-type conservation 
equations for both gaseous and droplet phases In two-dimensional, multiphase, 
combusting, turbulent flows. Nondlmenslonal 1 zatlon of these equations leads 
to the following two general forms (ref. 1 and private communication with X. 
Zhou and H.H. Chiu, University of Illinois at Chicago, 1983): 



where f x and f r are fluxes In the X and R directions, respectively; Tc 
is the transport coefficient for the variable E, Sc and S^ n t are the inner- 
phase and Interphase source terms, respectively, ana 1 takes the values 1=1 
for the cylindrical coordinate and 1=0 for the Cartesian coordinate. 

All conservation equations In these two forms are summarized In tables I 
and II, respectively. 

The physical models employed or developed In this study are Introduced 
In the following sections. 


Spray Spectrum Model 


The generalized Rosln-Rammler function is adopted to derive the size 
derivative of the droplet number density: 
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where S and t characterize the spray spectrum quality, r denotes gamma 
function, and n t and r^m are the droplet number density and the volume- 
mean Initial droplet radius. 


DROPLET Turbulent Diffusion Model 


Besides the trajectory motion, droplets also disperse through diffusion 
caused by the gas turbulence. A droplet diffusion model Is formulated (private 
communication with X. Zhou and H.H. Chiu, University of Illinois at Chicago, 
1983 and ref. 5) which takes the following final forms: 
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(5) 

( 6 ) 


(7) 


( 8 ) 


where og and a n i are the gas and droplet turbulent Schmidt numbers, 
respectively. 


Modified K-c Model 

The original K-c turbulence model Is modified to take account of 
droplet-phase effects. The droplets are assumed to share turbulent kinetic 
energy with the gas phase. The relationship of turbulent kinetic energies 
between gas and droplets can be approximated as 

~ K /* 2 

The modified K-equatlon Is expressed as 


v * (0p gV } " 7 * + eG k " ep g c ' V * J K A d *p (9) 

The e-equatlon remains unchanged. 
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Hybrid Turbulent Gas Combustion Model 
The combustion reaction Is supposed to take place In a single step: 


V + Y ox° 


Yp p 


( 10 ) 


A hybrid combustion model Is adopted; It assumes that the combustion rate 
Is controlled by the slower of two competitive rates of successive subprocesses 
the chemical reaction rate and the turbulent mixing rate. The former rate Is 
determined by the Arrhenius law, the latter rate Is calculated from the eddy- 
break-up ( EBU) model 

R p = - min ^ F>EBU I. I'V.Ar 1 ) (11) 


where 


R F,EBU = _C R ep g K 


1/2 



2 


1/2 


R F,Ar 


-Vg'gVox 


exp (-E/R u T g ) 


NUMERICAL METHODS 

A 42 by 38 staggered grid system Is superimposed on the flow domain. The 
previously mentioned two general forms of nonlinear partial differential equa- 
tions are then discretized Into their finite difference counterparts by using 
the control volume scheme and the primitive variables of pressure and velocity. 
The power law scheme (ref. 3) Is adopted to determine combined diffusion/ 
convection fluxes. The SIMPLER procedure is employed to predict gaseous flow 
fields, It reduces computer time by 30 to 50 percent In comparison with the 
popular SIMPLE procedure. The DROPLET procedure Is used for calculating drop- 
let flow fields. The alternating-direction double-reclprocatlng-sweep line- 
by-line TDMA method Is adopted to solve simultaneous algebraic equations Iter- 
atively. The block correction technique Is used to facilitate convergence. 

The usual wall function method Is employed to bridge the near wall regions 
where the laminar viscosity effect Is quite strong. This method Is very effec- 
tive In greatly reducing computer time while obtaining rather satisfactory 
results. All these methods have been assembled Into the computer code GEMCHIP. 
The main flow chart of GEMCHIP Is shown In figure 1. 

A droplet-free solution Is first obtained through SIMPLER, which provides 
both final results for the single gas-phase flow and the Initial guesses of the 
gas flow fields for multiphase flow. Then we enter DROPLET to calculate drop- 
let flow fields and the Interphase exchanges of mass, momentum, and energy. 
These quantities are substituted as source terms Into SIMPLER, and the gas flow 
fields are modified. Then we enter DROPLET again. This alternating Iterative 
process continues until all the gas- and droplet-flow equations satisfy the 
convergence criterion In all grids. The typical Iteration number ranges from 
40 to 75 for gas-phase burning problems and from 120 to 150 for multiphase 
burning problems. The corresponding computer times are about 5 and 30 min, 
respectively, on an IBM 4342 machine. 
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THE AFTERBURNER STUDIED 


Figure 2 shows a low-bypass-ratio turbofan afterburner In which the com- 
puter modeling Is made. The afterburner consists of a cone-type diffuser/mixer, 
a cylindrical combustion section, and a single-ring bluff-body flameholder. A 
cylindrical partition board separates the cold fan flow from the hot core flow. 
The temperature and velocity are 350 K and 80 m/sec In the fan flow and 600 to 
800 K and 60 to 200 m/sec In the core flow. The liquid fuel Is evenly Injected 
In spray form Into the core flow, and then the droplets are heated, vaporized, 
mixed, and combusted with oxygen In the airflow. 


Handling of Irregularly Shaped Wall and Flameholder 

The curved afterburner wall and the flameholder are simulated by blocking 
off some of the control volumes of the regular grid system so that the remain- 
ing active control volumes form the desired shape. The advantage of this 
method Is that a computer code developed for the nominally regular domain can 
be adapted to flows of arbitrarily shaped configurations so that there Is no 
need to generate new programs and the time for developing and testing such 
programs Is saved. 

The Irregular boundary can be represented approximately by a series of 
rectangular steps (fig. 3). Usually a rather crude representation of the 
boundary can produce surprisingly good results. 

The values of the relevant flow parameters should be given In the blocked- 
off regions. This Is done by using extremely large source terms or by assign- 
ing extremely large transport coefficients In these Inactive control volumes. 
These Inactive grids also participate In the Iteration process. 


Boundary Conditions 

Appropriate boundary conditions are selected to represent the wall regions 
of the afterburner. The following assumptions can be realistically made with 
regard to this kind of selection: 

(1 ) No slip at the wall 

(2) Zero turbulent fluctuation at the wall 

(3) Adiabatic and Impermeable wall 

It Is reported that the predicted solution Is very sensitive to the Inlet 
conditions. For this reason, the Inlet profiles should conform to measured 
values whenever available. Unfortunately, measurements of the profiles of the 
Inlet variables, especially the Inlet velocity vectors and the Inlet turbulence 
variables K and e, are seldom carried out In complex geometries. It Is, 
therefore, reasonable to assume the following profiles (ref. 6): 

U „ = U (1 - r) n 
g go' 

v = o 
g 
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K - K 0 (l ♦ m(r) 2 ) 


* - * 0 (^ o ) 1 ' 5 

where 0.2 > n > 0 and 2 > m > 0. The variables with subscript o denote 
values at the centers of core flow or fan flow. The nondlmenslonal radius 
r = r/R for the core flow and r = |r - (R-| + R 0 )/2|/((R-| + R 0 )/2) for the fan 
flow, where R 0 denotes the radius of the partition board and Ri represents 
the diffuser's Inlet radius. 

The Inlet profiles of the droplet number densities, temperatures, and 
velocities vary with the fuel nozzle type and spray quality. 

In the axis of symmetry the radial gradients of all droplet-phase vari- 
ables are taken to be equal to zero. 

The exit plane should be located far downstream and outside the recircula- 
tion zone so that the local parabolic flow assumption applies and that the 
calculating domain Is Isolated from the ambient environment. The only excep- 
tions are the equations for pressure p and for the pressure correction param- 
eter p', since the pressure effect Is always two way. The problem Is solved 
by calculating the outlet velocity component In the main flow direction. 


RESULTS AND DISCUSSIONS 

The objectives of the numerical study are to analyze the shape and the 
size of recirculation zones at different working conditions, to examine the 
vaporization process of multisize fuel sprays, to compare the multiphase flame 
structures with the gas flame structures, and to study the parametric sensitiv- 
ity of afterburning efficiency. 


Recirculation Zones 


The gas velocity vector fields Inside the afterburner are shown In 
figures 4 to 6 for gas flow and multiphase flow. The numerical calculations 
reveal that the recirculation zones In the wake of the flameholder are Induced 
by negative pressure gradients, which. In turn, are caused by gaseous viscosity. 


Figures 4 and 5 are pictures of velocity vectors and X-dlrectlon velocity 
component profiles In a gas phase flow. The length of the recirculation zone 
In the cold flow case Is about 2.0 to 3.5 times the flameholder width. It Is 
readily seen from figure 4 and 5 that the wake effect still exists far down- 
stream of the flameholder. Some experimental reports claim that the reclrcula- 
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expansion and the reduction of the absolute value of the negative pressure 
behind the flameholder. But this trend Is not obvious In this calculation. 
Research work Is under way In this direction. 


Figure 7 shows droplet velocity vectors of two size groups In the multi- 
phase afterburning condition. A comparison of the multiphase flow with the 
gas flow field (fig. 6) shows that the small droplets are capable of reaching 
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velocity equilibrium more quickly than the large droplets. The velocity non- 
equilibrium In the downstream Is somewhat enhanced because the gas Is acceler- 
ated more rapidly than the droplets. This Is especially true for the large 
droplets. The droplet recirculation zone Is obviously smaller than the gas 
recirculation zone. 


Cascading Vaporization 

The preheating and vaporization processes of droplets of different size 
groups are Illustrated In figure 8. It Is readily seen from the predicted 
curves that the small droplets (e.g., D-j = 0.15) complete the preheating 
process and Initiate the vaporization process much earlier than the large 
droplets (e.g., 0] = 1.60). 

The phenomenon of successive Initiation of the vaporization of droplets 
of different sizes In spray Is designated as cascading vaporization (ref. 1) 
since the Initial vaporization lines of different size droplets constitute a 
form of cascade. 


Gas Flame and Multiphase Flame Structures 

Figures 9(a) and (b) show the temperature field and flame structure of 
the gas combustion flow. The gas flame structure Is rather simple: the flow 

of the combustible mixture mixes with the fan flow of the cold air, and the 
flame Is stabilized In the recirculation zone close behind the flameholder. 
There Is a "dark" area, which Is fuel-deficient, In the wake of the flame- 
holder. This Is where the combustion Is completed. 

The structure of the multiphase flame Is very complicated. Calculation 
reveals that there are three principal combustion modes: (1) pre-evaporating 

flame, (2) dense spray group flame, and (3) dilute spray group flame. The 
total fuel-air ratio and the spray group combustion number G c are the two 
main parameters determining combustion mode. The G c Is actually the ratio 
of the characteristic vaporization rate to the convection flow rate and was 
first proposed by Prof. H.H. Chiu (ref. 7). 

Pre-evaporating flame . - The pre-evaporation flame Is shown In figure 10. 
When the G c number Is very large (fine atomization, l.e., the spray con- 
sists of an extremely large quantity of small droplets) all the droplets are 
vaporized at a typical afterburner Inlet temperature before they get to the 
flame zone. The flame structure Is similar to that of gas phase combustion. 

Dense spray group flame . - If the G c number Is rather large and the 
fuel-air ratio Is also high, dense spray group combustion occurs (fig. 11); 
this Is characterized by the presence of an oxygen-deficient zone, which 
decreases combustion efficiency. Both types of gas and droplet burning coexist 
In the flame. 

Dilute sprav group flame . - Thlsjnode occurs at a small G c number, which 
Is characterized by poor atomization quality and a low fuel-air ratio: that 

Is, there Is a small quantity of big droplets. If the droplet size Is large 
enough and droplets concentrate near the central line of the afterburner, the 
v-type flame shown In figure 12 may appear. This type of flame Is unstable. 
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Figures 13(a) and (b) show, respectively, the axial distributions of 
afterburning efficiencies In gas-phase and multiphase combustion cases. It Is 
seen that the growth of the combustion efficiency for the multiphase flame has 
a slower start than that for the gas flame. This difference Is caused by the 
droplet preheating and by vaporization processes. The results also Indicate 
that the smaller the spray group combustion number G c , the larger the 
difference. 


Parametric Sensitivity Study of Combustion Efficiency 

Figure 14 shows the effect of fuel-air ratio on gas-phase combustion effi- 
ciency at two different pressures (2.0 kg/cm 2 and 0.5 kg/cm 2 ). The effi- 
ciency curve reaches Its peak at an approximately stoichiometric fuel-air ratio 
In the core flow. At higher pressures, the efficiency peak moves toward richer 
mixtures and the curve becomes flatter because of the dilution action by the 
fan flow and the Improvement of combustion conditions. 

Figure 15 shows the variation of the multiphase combustion efficiency with 
the fuel-air ratio at two different spray group combustion numbers. With 
multiphase combustion there Is a combustion efficiency peak at an appropriate 
fuel-air ratio similar to that with gas-phase combustion. With an Increase of 
G c number, the efficiency peak moves toward leaner mixtures - which means 
mixtures with a larger number of smaller-sized droplets will have a higher com- 
bustion efficiency at a lower fuel-air ratio. The results are also In agree- 
ment with the experimental measurements. 

Figures 16(a) and (b) show the effect of the Inlet gas temperature on com- 
bustion efficiencies of gas-phase and multiphase flames. The temperature 
Increase enhances the heat release rate and speeds up droplet vaporization, 
which. In turn, enhances the combustion efficiencies of both the gas and multi- 
phase flames. But, for obvious reasons, the temperature effect Is stronger In 
multiphase flames. 

The effect of Inlet gas pressure on combustion efficiency Is similar to 
that of Inlet temperature. It Is seen from figure 14 that. In the case of gas 
phase combustion, an Increase of gas pressure from 0.5 to 2.0 kg/cm 2 signifi- 
cantly Improves combustion. Figure 17 shows the Increase of combustion 
efficiency with gas pressure In a multiphase combustion case. The pressure 
effect Is negligible at pressures greater than 1.3 kg/cm 2 . 

Figure 18 shows the effect of gas-flow velocity on gas-phase combustion 
efficiency. An Increase In gas velocity results In a decrease In combustion 
efficiency because of the reduction of the residence time of the gas flow In 
the afterburner. The higher the velocity, the lower the efficiency. The 
effect of Inlet gas velocity on multiphase combustion efficiency Is rather 
complicated since the velocity affects not only the residence time but also 
the droplet atomization, vaporization, and spray group combustion number G c . 


CONCLUSIONS 

The mathematical models and numerical methods for predicting multiphase 
turbulent reacting flows, developed In the authors' study of spray group com- 
bustion phenomena have been successfully extended and applied to the numerical 
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study of flow fields and combustion characteristics of a low-bypass-ratio 
turbofan afterburner. The computer code developed by the authors Is versatile 
and effective In this study. The technique of blocking off relevant control 
volumes makes It possible for a computer code developed for a regular grid 
system to be applied to arbitrary flow domains. 

The numerical analysis reveals that the size of the recirculation zone 
behind the flameholder Is about 2.0 to 3.5 times Its width. The droplet recir- 
culation zone Is smaller than the gas recirculation zone and decreases or even 
disappears with an Increase In droplet size. 

The calculation also exposes the existence of three principal combustion 
modes In multiphase combustion. They are the pre-evaporating flame mode, the 
dense spray group flame mode, and the dilute spray group flame mode. Study of 
flame structure Is useful to afterburner design optimization. 

The results reveal the dependence of afterburner combustion efficiency on 
the fuel-air ratio and on the operating parameters (such as the Inlet gas tem- 
perature, pressure, and velocity). As for multiphase combustion, which 
prevails In turbofan afterburners, the spray group combustion number has a 
significant effect on combustion efficiency and flame structure. It Is con- 
cluded that afterburners should be designed to work In the optimum spray group 
combustion number, which Is related to the fuel-air ratio and Its distribution. 
Generally speaking, the favorable design has a pre-evaporating flame mode or a 
dense group combustion mode with an Initial fuel distribution capable of keep- 
ing oxygen-deficient zones as small as possible. 

Although the previously predicted results are In good qualitative agree- 
ment with the known facts, accurate diagnostic measurements and detailed 
experimental verifications are necessary to examine the quantitative agreement 
and to Improve the present models. 
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TABLE I. - THE GOVERNING EQUATIONS IN THE FORM OF EQUATION (1) 
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TABLE I. - Continued. 
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*G-equat1on used for turbulent diffusion flame; not mentioned In text. 




















(b) Form of equation (2) 



go go 




















Figure 1. - Main flowchart of GEMCHIP. 
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Figure 3. - Blocked-off region In regular grid system. 



Figure 4. - Gas velocity vectors of a cold flow. 
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(a) Temperature field In a gas combustion flow. 



(b) Gas fuel concentration field In a gas phase combustion flow. 


Figure 9. - Gas phase combustion flow. 
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Figure 10. - Pre-evaporating flame mode. Droplet preheating and 
evaporating zone, (1); gas phase mixing zone, (2); gas phase burning 
zone, (3); fuel-deficient zone, (4). 



Figure 11. - Dense spray group flame mode. Droplet preheating and 
evaporating zone, (1); gas phase mixing zone, (2); gas phase burning 
zone, (3); gas phase and droplet burning zone, (4); oxygen-deficient 
zone, (5). 



Figure 12. - Dilute spray group flame mode. Droplet preheating and 
evaporating zone, (1); gas phase and droplet burning zone, (2); gas 
phase combustion zone, (3); fuel-deficient zone, (4); gas phase 
mixing zone, (5) . 
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Figure 14. - Relation of gas phase combustion efficiency 
with fuel-air ratio and pressure. 
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Figure 15. - Relation of multiphase combustion efficiency 
with fuel-air ratio and spray group combustion number. 
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Figure 16. - Effect of Inlet temperature on 
combustion efficiency. Gas phase flame, 

(A) ; multiphase flame (at better condition), 
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Figure 17. - Effect of inlet pressure on multiphase 
combustion efficiency. 
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Figure 18. - Effect of Inlet gas flow velocity on gas 
combustion efficiency. 
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MODELING TURBULENT, REACTING FLOW 

Russell W. Claus 
NASA Lewis Research Center 
Cleveland, Ohio 


This paper examines several of the approximations or models Involved In 
the development of a numerical combustor flow code. In the first section, the 
Importance of numerical accuracy Is Illustrated, and the Impact that Improved- 
accuracy schemes have on slowing convergence Is demonstrated. Solution algo- 
rithms that can speed convergence are discussed and some performance features 
of these algorithms are Illustrated. A sample calculation displaying the 
Importance of boundary conditions on a three-dimensional numerical prediction 
Is presented. The Inaccuracy of a current turbulence model In highly turbu- 
lent (nonequilibrium) regions Is described. Finally, the surprisingly good 
performance of a six-flux model In describing radiation heat transfer Is dis- 
played. In all the areas examined, continued research Is still needed, but 
valuable engineering tools are available today. 

INTRODUCTION 

Three-dimensional combustor calculations Involve detailed modeling of 
several Important physical processes. Airflow, chemical reactions, fuel 
sprays, and turbulence are just a few of the physical processes that must be 
described. Many of these processes occur on both a molecular and a macro- 
scopic scale. To exactly describe these processes numerically, one must 
resolve these scales on a computational mesh. And this Is clearly beyond 
current computational resources. To make the computational task tractable, 
we Introduced modeling assumptions. These modeling assumptions limit the 
generality of the computational flow code, but It Is hoped that the dominant 
physics remain correctly represented. 

Modeling assumptions are only the first limit of generality Introduced 
when developing a combustor flow code. A further limit Is Introduced by the 
need to approximate the modeled equations before they are solved numerically. 
This approximation process can significantly affect the accuracy of a model 
prediction. 

In view of all these factors, a combustor designer cannot be expected to 
fully embrace a computer model prediction. If a clear distinction between 
modeling errors and errors Introduced by the numerical solution algorithm can- 
not be made, then the designer Is left with a very unreliable computational 
tool . 


This paper examines some of the compromises made when developing combus- 
tor flow codes and how these compromises affect the accuracy of turbulent flow 
calculations. To assess the balance of the modeling compromises, we compare a 
number of sample calculations with experimental data. This paper examines the 
five areas of numerical accuracy, solution algorithms, Inlet boundary condi- 
tions, turbulence models, and radiative heat transfer. 
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The NUMERICAL ACCURACY section focuses on the development of more accu- 
rate numerical methods to be used In combustor flow codes. Upwind differenc- 
ing, which Is currently used in these flow codes, Introduces an appreciable 
error (or numerical diffusion) Into the calculation. This error may be of 
such a large magnitude that It obscures the turbulence model used In the cal- 
culation. A series of calculations are Illustrated which demonstrate the 
accuracy of a variety of differencing schemes. 

An Important aspect of Improved accuracy Is the effect that the differ- 
encing scheme has on the rate of convergence. The Improved accuracy schemes 
all appear to require more CPU time to converge. For this reason, the next 
section discusses solution algorithms. SIMPLE ( seml-lmpl 1 c 1 t pressure-linked 
equations) Is one of the most widely used solution algorithms for solving the 
steady-state form of the Navler-Stokes equations. Although this scheme has 
proven to be quite effective. Its convergence rate can be Improved. This 
paper focuses on two approaches which accelerate convergence by performing 
corrections that Improve the Iterative agreement with continuity. These 
alternate schemes are Illustrated In a series of calculations. 

A third section examines the Importance of Inlet boundary conditions. An 
Illustrative example Is displayed. 

The fourth section discusses turbulence models and reaction closures. 

Here the more pragmatic approaches to calculating turbulent reactions are 
Illustrated. An eddy-breakup model and a PDF (probability density function) 
method are described and compared. 

Finally, the fifth section examines the accuracy of a radiation heat- 
transfer model. A six-flux model of radiative heat transfer Is described. 

This model provides only a limited geometric description of the radiation 
transfer process, but a comparison with experimental data Indicates an encour- 
aging level of agreement. 


NUMERICAL ACCURACY 

The following Is a general form of the time-averaged equations that must 
be solved In a combustor flow code: 

f^(pU<p) + f^(pV<p) - fj) - fj) = (1) 

where <p can represent U, V, W, u^Uj, K, e, or H; where Is the 

diffusion coefficients (e.g., y e ff); and s <p the source term (e.g., - 3P/3x). 
The equation represents convection of a conserved scalar that Is subtracted by 
diffusion terms equal to a source term that could represent a pressure gradi- 
ent or a source/sink term. These equations are solved by discretizing on a 
staggered mesh using the finite volume method (fig. 1). The staggered mesh 
system Is used to avoid the pressure-velocity decoupling that can result In a 
finite volume representation of Incompressible flow. 

An example of the discretization of a convective term using central dif- 
ferencing Is the following: 
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— — f f dx dv 

AxAy J J dx y 


( 2 ) 




(pu«p) e - (pu«p) w 
AX 


(Alternate types of differencing are possible and some will be Illustrated In 
the following paragraphs.) Once all the terms have been approximated, like 
terms are then rearranged In a substitution formula which can be solved using 
a tri-diagonal matrix algorithm (TOMA) procedure applied In alternate 
directions (ADI): 

*p<Pp « a N «p N + a s «p s ♦ a E q> E + <ypg ♦ S (3) 

where, for example. 


a N = 




The source and diffusive terms In these equations are approximated using 
central differencing which Is second-order accurate (Error « 0(ax)2). The 
convective terms have typically been differenced using hybrid differencing 
which reverts to upwind differencing when the absolute value of the cell 
Reynolds number, or cell Peclet number. Is greater than 2, which Is the main 
reason for the loss of numerical accuracy. An example of convective discreti- 
zation using hybrid differencing Is as follows: 


au<p 

dx 




(U<p) - ( u<p) 

e w 

Ax 


(4) 


For upwind, with u > 0, 


( u«p) = u <p_ 
e e P 


and 


<“*>„ - 2 u w<*w * V 

Since upwind differencing Is only first-order accurate. It can Introduce an 
extensive amount of numerical diffusion Into the calculation. The virtue of 
this scheme Is that It provides "bounded," stable solutions. Higher order 
convective schemes Invariably yield nonphysical oscillations In the solution. 
Upwind differencing avoids this, at a considerable loss In accuracy. 

To alleviate this problem, NASA has conducted a program to Identify and 
Incorporate an Improved accuracy differencing scheme Into a combustor flow 
code. Under a portion of this program a variety of differencing schemes were 
examined In several test calculations. The schemes examined Included QUICK 
(quadratic upstream Interpolation) and SUD (skewed upwind differencing). 

QUICK differencing was developed by Leonard (ref. 1). This scheme Improves 
the accuracy of convective differencing by performing an upwind biased quad- 
ratic Interpolation. For u > 0, 
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( 5 ) 


(U») e . u e (- 3 » w * f„p . §* E ) 

<“*>« - u w (- e^ww * «*w * 8»p) 

where grid point locations are as noted on figure 2. This scheme Is second- 
order accurate and can produce nonphysical oscillations In the solution. SUO 
(skewed upwind differencing) (ref. 2) attains high accuracy by differencing In 
an upwind manner along the flow streamlines. While maintaining the same for- 
mal accuracy as upwind differencing, the truncation error In SUD Is smaller. 

For example: 

For u > 0 and V > 0, 

(u<p) e = VP | 

(u <p) w = u w O " a)<p w + u w a<p sw ( (6) 

V = minimum of (l,V/2u) I 

where grid point locations are as noted In figure 3. As with QUICK, SUD can 
produce nonphysical oscillations In the solution; therefore, a scheme to 
"bound" SUD was also examined. This scheme employs the concept of flux- 
blending (ref. 3), wherein a bounded flux determined from upwind differencing 
Is blended with the unbounded, but more accurate, SUD flux. The main factor 
Is to blend as little of the lesser accurate scheme while still maintaining a 
properly "bounded" solution. This procedure, called BSUDS, starts from an 
Initial, totally skew-differenced estimate and blends an upwind flux If the 
solution Is out of the range of neighboring values. If the solution Is In 
range (l.e., bounded), then no blending Is performed. 

An Illustration of the accuracy of the upwind, QUICK, and SUD schemes Is 
seen In figure 4. This figure displays the results of a single-point scalar- 
transport calculation made for various flow angles. All schemes agree with an 
exact solution (no error) at a zero flow angle; however, at angles greater 
than zero, each scheme displays some degree of error relating to numerical 
diffusion. The error displayed by upwind differencing Increases with flow 
angle to a maximum at 45°. The QUICK scheme displays a similar behavior, but 
the overall error Is much less. The SUD scheme displays a maximum error 
around 15°, but It tends to zero at angles approaching 45°. Both QUICK and 
SUD display a much higher level of accuracy than upwind. 

Although a scalar transport calculation Is useful for a general examina- 
tion of some aspects of differencing scheme performance, a laminar flow calcu- 
lation Is a more complete test. The results of a series of laminar flow 
calculations from reference 4 are displayed In figure 5. In this figure, 
axial velocity profiles at a distance of one-half a duct height from the Inlet 
are shown for two different computational meshes. In these calculations, the 
steepness of the velocity profile Indicates accuracy. Steep velocity profiles 
are exhibited by QUICK and BSUDS; the upwind profiles exhibit a high degree of 
numerical diffusion. On the coarse mesh, BSUDS appears to be more accurate 
than QUICK; on the fine mesh, there Is not much of a distinction between the 
two schemes. 
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An Important aspect of Improved accuracy Is the effect that the differ- 
encing scheme has on the rate of convergence. The computational times to 
converge the system of governing equations In the previous laminar calcula- 
tions are shown In table I. The convergence times are ratloed to the upwind 
convergence times to clearly Illustrate the computational penalty paid to 
attain Improved accuracy. Generally, the Improved accuracy schemes required 
from 3 to 15 times longer to reach a converged solution. To a degree one can 
hope that this computational penalty can be offset by using coarse meshes to 
achieve the same overall level of accuracy. In reality, a relatively fine 
mesh Is needed, even with the high-accuracy schemes. In any case, the need 
for Improved solution algorithms for these more accurate differencing schemes 
Is strongly Indicated. 


SOLUTION ALGORITHMS 

The previous section demonstrated the need for Improved solution algo- 
rithms. This section examines the widely used SIMPLE algorithm and two modi- 
fications to this scheme. These schemes cover only a small portion of the 
wide range of methods to accelerate solution convergence. Vectorlzatlon, 
direct-solution methods, and multigrid methods - to list just a few - are all 
areas of active research that are certain to yield much greater computational 
benefits In the near future. 

In the SIMPLE algorithm, a guessed pressure field Is Inserted Into the 
discretized momentum equations to obtain a velocity field. The pressure field 
Is corrected by an equation which Is derived through a combination of the con- 
tinuity and momentum equations. The velocity field Is then updated and Is 
used In the solution of the equations for k, e, and <p. The corrected pres- 
sure field Is treated as the guessed pressure field, and the procedure Is 
repeated until a converged solution Is obtained. 

The following velocity correction equation Is used In the SIMPLE scheme 
for u at point e: 

A u' =5" A . U' + ( P * - P ' ) A 
e e ^ nb nb P E e (7, 

Primes Indicate corrections to old values. The underlined term, which repre- 
sents the Influence of corrected pressures on neighboring velocities. Is 
neglected In the SIMPLE algorithm. The converged solution Is unchanged by the 
exclusion of this term since, for the steady-state solution, the corrections 
go to zero. Neglecting terms, however, does force the use of low underrelaxa- 
tlon factors, which can slow convergence. 

The SIMPLER algorithm Improves on the SIMPLE scheme by Including the 
previously neglected terms when calculating the pressure field. The calcula- 
tion sequence starts with a guessed velocity field. An equation that solves 
for the pressure field (using the terms Ignored In the SIMPLE scheme) Is cal- 
culated from the guessed velocity field. This pressure field Is then used to 
solve the discretized momentum equations to obtain a velocity field. The 
velocity field Is corrected In a manner similar to the SIMPLE velocity correc- 
tion. This velocity field Is then treated as the guessed velocity field, and 
the Iteration procedure Is repeated until convergence Is reached. 
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Because additional equations are solved, each Iteration through the 
SIMPLER routine Involves more computational time than an Iteration step through 
SIMPLE. However, higher underrelaxatlon factors can be applied In the SIMPLER 
routine, thereby accelerating convergence. 

The PISO scheme also takes Into account the terms neglected In the SIMPLE 
code, but In a different manner. The PISO routine mimics the SIMPLE approach 
until the end of the first Iteration. At this point, the PISO scheme employs 
an equation containing the neglected terms to correct the pressure and veloc- 
ity field to more closely agree with continuity. Again, this procedure Is 
repeated until the solution converges. In this manner, the PISO code allows 
for higher underrelaxatlon factors to accelerate convergence. 

The performance of these various solution schemes Is displayed In 
figure 6. As an example, the convergence times for a 38-by-38 grid point 
calculation are plotted as a function of underrelaxatlon factor (fig. 6(a)). 
When a reasonably large underrelaxatlon factor Is used, SIMPLER or PISO con- 
verge about twice as fast as SIMPLE. In addition, SIMPLER and PISO converge 
over a larger range of underrelaxatlon factors than SIMPLE. From an engineer- 
ing standpoint, Improving the "robustness" of a computational scheme Is often 
just as Important as accelerating convergence. The computational benefit of 
using SIMPLER and PISO for fine-mesh calculations Is displayed In figure 6(b). 
The greater the number of mesh points used In the calculation, the greater the 
benefit of SIMPLER or PISO over SIMPLE. For example, a calculation of approx- 
imately 3300 grid points converges three times faster using PISO or SIMPLER 
than It does using SIMPLE. This Is a savings of about 600 CPU seconds. The 
1440 grid points calculation demonstrated a savings of only about 100 CPU 
seconds . 


INLET BOUNDARY CONDITIONS 

In any calculation of a complex, three-dimensional turbulent flow, the 
boundary conditions that are needed In the calculation are frequently unknown 
or unmeasured. The use of Inappropriate values at the computational boundary 
can sometimes be the main limit to the calculation's predictive capability. 
This error can sometimes be more Important than numerical accuracy or turbu- 
lence model considerations. Figure 7 shows an example of a three-dimensional 
jet-ln-crossf low calculation using alternate boundary conditions. In one 
calculation, the jet orifice flow was specified as having a uniform plug flow 
at a position two jet diameters upstream of the orifice outlet. This allowed 
the flow to distort as It exited the orifice outlet. The second calculation 
specified a uniform plug flow at the orifice outlet. The resulting axial 
velocity profiles are compared with experimental data In figure 8. Somewhat 
surprisingly, the uniform boundary condition at the orifice compared more 
favorably with experimental data than the theoretically more correct distorted 
profile. This should not be Interpreted as an endorsement of the use of 
uniform-plug-flow boundary conditions for these types of flows. Indeed, there 
Is some Indication of experimental error (ref. 5). This example Is meant only 
to Illustrate that unmeasured or unknown boundary conditions can significantly 
affect a numerical calculation. 


TURBULENCE MODELS 

Although a great deal of progress Is being made In solving the three- 
dimensional, time-dependent Navler-Stokes equations In large-eddy or direct 
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numerical simulations, practical engineering calculations currently require 
the introduction of some form of turbulence modeling. These models are based 
on either Reynolds or Favre averaging of the exact Navier-Stokes equations, 
reducing the unsteady form of these equations to an averaged form. Currently, 
the most widely used turbulence model is the two equation, k-e, closure. This 
model relates the Reynolds stresses to a turbulent viscosity through 
Boussinesq's eddy-viscosity concept: 


-pu 1 u j = pv 


faUj 
T\ixj 
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where vj = C u k 2 /e. The turbulent viscosity is related to the kinetic 

energy k and the dissipation rate c of the turbulence. The transport 
equations are then solved for k and e: 
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Where C = 0.09, C _ = 1.44, C _ = 1.92, a. = 1.0, o = 1.3, a = 0.9, k 

v el cl K e <p 
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The model constants typically employed are those recommended in 
reference 6. The turbulent Schmidt number a ip , is frequently changed from 
0.9 to as low as 0.2, depending on the flow being studied. 


This two-equation model is based on several assumptions which should be 
considered when making a numerical calculation. First, the flow is assumed to 
be close to equilibrium; that is, the flow properties change relatively slowly. 
Second, the turbulence Reynolds number is assumed to be high. Third, the tur- 
bulence is assumed to be isotropic. 


The main concern is how well this model, with its Inherent assumptions, 
can represent combustorlike flow fields. Figure 9 displays a comparison 
between laboratory experiments and numerical predictions of two different iso- 
thermal flows. (Figs. 10 and 11 show the locations of the measurements that 
were made in the flow fields.) In the two-dimensional bluff-body comparison, 
a major disagreement between measurements and predictions is evident at an 
axial distance x/D of approximately 0.8. This corresponds to the end of the 
recirculation zone and causes an Incorrect prediction of recirculation zone 
length. The jet-in-crossf low comparison displays a similar disagreement in 
the region where the turbulence intensity is high. In this comparison, both 
hybrid and BSUDS differencing were used in the predictions. It is obvious 
from the displayed results that numerical accuracy can have a major Impact on 
the comparison with experimental measurements. Hybrid differencing is so 
completely Influenced by numerical diffusion that the qualitative agreement 
between experiment and calculation, evident in the BSUDS results, is elimi- 
nated. The BSUDS results are not grid-independent, but it seems unlikely that 
this will fully explain the noted disparity. 
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Both of these flow fields display the greatest disagreement between exper- 
iment and calculation where the turbulence Intensities are the highest. These 
are regions where the flow field Is likely to be far from equilibrium. It Is 
Interesting to note that a full Reynolds stress transport (RST) model calcula- 
tion (presented by McGuIrk, J.3., Papadlmltrlou, C., and Taylor, A.M.K.P. at 
the Fifth Symposium on Turbulent Shear Flows held at Cornell University, 

Ithaca, New York, August, 1985) did not yield appreciably better results for a 
similar flow field. Because both models appear to lose validity around the 
region of the stagnation point, further model development Is needed. 

Reaction closures Involve a further series of assumptions and approxima- 
tions. The simplest level of closure Is to assume that the reaction Is mixing 
limited and Ignore the effect of temporal density fluctuations. The reaction 
rate can then be related to an eddy lifetime (e/k) using either a Magnussen- 
Hjertager or Spaldings eddy-breakup model: 

Rate = -S (10) 

m p 

where 

T Mo 2 M r 1 

S m p = " p k m1n | Am F» A STOIC ’ AB STOIC + 1J 


or 


S m p “ " C R p k m Fu 

where A, B, and Cp are empirical constants and e/k can be considered as 
the eddy lifetime. The density and temperature throughout the flow field can 
then be established from equilibrium chemistry. 

Improved physical realism can be added to this model by Introducing the 
probability density function (PDF) for mixture fraction f. This can account 
for the unsteady time history of the mixture fraction at each point In space. 

A nonlinear functional dependence between concentration, temperature, and 
density occurs as pockets of alternately hot and cold gases are swept past 
each point In space. Integrating the resultant PDFs allows a determination of 
the mean properties of the flow. 

Figure 12 displays the results of both an eddy-breakup and a PDF model 
calculation of a planar-mixing layer. The PDF model calculations (ref. 7) 
demonstrate a significant Improvement over the eddy-breakup results. The main 
factor Involved In this Improvement Is the more physically correct representa- 
tion of the mixture fraction variation at the point of maximum temperature 
rise. Of course, both of these models are only as good as the turbulence 
closure and are wholly Inadequate when finite-rate chemistry Is Important. 

RADIATION HEAT TRANSFER 

The final topic to be covered In this review concerns radiative heat 
transfer. Radiation Is one of the most significant and least understood heat 
loads to the combustor liner. It can also play a significant role In the 
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determination of flame temperature. Any numerical description of a gas tur- 
bine combustor must Include a radiation heat-transfer model. The six-flux 
model Is the model most commonly used to approximate multidimensional radi- 
ative transfer. In this model, differential equations describing the radi- 
ative fluxes In positive and negative directions along the principal axis are 
solved: 


d ( 1 dR x \ 

dx\a + S dx / 

1 d / r dR r 

r dr \a ♦ S ♦ 1 dr 

1 d ( 1 dR Z \ 
r de\a + S r de J 


a(R x - 
a(R 

= a(R z 


E) + |(2R X - R r - R Z ) 
r - E) + |(2R r - R X - R Z ) 

- E) + |(2R Z - R X - R r ) 


(ID 


where 

RX,r,z composite fluxes 
a absorption coefficient 

S scattering coefficient 

E oT 4 

The main Input to this analysis concerns the optical characteristics of the 
hot gas and soot which must be arbitrarily specified or calculated through a 
soot formation and oxidation model. 

The performance of the six-flux model (ref. 8) Is displayed In figure 13. 
Although the model overestimates the radiative heat transfer In comparison 
with experimental data, the qualitative trend Is quite closely followed. Given 
the large number of approximations used In the analysis, the agreement with 
experimental data Is quite surprising. The six-flux model does not accurately 
treat the angular dependence of energy transfer, and the determination of the 
optical characteristics of the soot cloud still remains as an area of needed 
research; however, fairly good results appear to be possible In this example. 

CONCLUDING REMARKS 

Perhaps the most Important question any review on numerical modeling can 
answer Is whether or not current computational codes can be usefully employed 
In the design of combustion devices. Certainly a great deal of research Is 
needed before one can expect quantitative predictive accuracy, and It seems 
likely that some hardware problems will only be resolved through development 
testing. The best computer program will never replace the designer's Innova- 
tive mind, but computer predictions can be used to extend the designer's pro- 
ductivity. New designs can be examined much more rapidly on the computer than 
In hardware testing. Development costs can be reduced. The promise of this 
computer-based design methodology Is so great that these numerical models will 
be used despite their deficiencies. Designers should not and probably will 
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not abandon empirical design tools, but the cautious adoption of numerical 
models In the design process Is a trend which can reap Important benefits. 
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TABLE I. - RATIO OF CONVERGENCE TIMES 
FOR VARIOUS DIFFERENCING SCHEMES 
WITH UPWIND CONVERGENCE TIMES 
USED AS THE STANDARD 


Mesh 

Upwind 

BSUDS 

QUICK 

Coarse (30 by 22) 

1 

6.4 

A 

0 . c 

Fine (58 by 38) 

1 

14.7 

15.7 
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Figure 1. - Staggered mesh system for discretizing equations 
by finite volume method. 
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Figure 2. - Staggered mesh system for 
quadratic upstream interpolation. 



Figure 3. - Staggered mesh sys- 
tem for bounded skew upwind 
differencing scheme (BSUDS). 
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Figure 6. - Convergence times. 
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(b) For various computational meshes. 
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Figure 7. - Three-dimensional jet-in-crossflow 
calculation. 
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Figure 8, - Calculated axial velocity profiles com- 
pared with experimental data. 
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(a) Two-dimensional bluff-body wake. (b) Three-dimensional jets in crossflow. 

Figure 9. - Comparison of laboratory experiments with numerical predictions of two different isothermal flows. 
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Figure 10. - Two-dimensional bluff-body geometry. 
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(b) Eddy breakup model. 


(c) PDF model. 


Figure 12. - Results of calculations of a planar-mixing layer. 
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(a) Schematic drawing of experimental measurement 
locations. 



(b) Performance of six-flux 
model. 

Figure 13. - Radiative heat transfer test calculation. 
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EXPERIMENTAL AND ANALYTICAL INVESTIGATION OF THE VARIATION OF SPRAY 

CHARACTERISTICS ALONG A RADIAL DISTANCE DOWNSTREAM OF A 

PRESSURE-SWIRL ATOMIZER 

J.S. Chin, W.M. LI, and X.F. Wang 
Beijing Institute of Aeronautics and Astronautics 
Beijing, People's Republic of China 


The variation of spray characteristics along a radial distance downstream 
of a pressure-swirl atomizer was measured by laser light-scattering technology. 
An analytical model was developed to predict the variation of spray character- 
istics along the radial distance. A comparison of the predicted and experi- 
mental data showed excellent agreement. Therefore, the spray model proposed, 
although relatively simple. Is correct and can be used, with some expansion 
and modification of the prepared model, to predict more complicated spray 
systems . 


INTRODUCTION 

The characteristics of the spray formed by an atomizer In a liquid-fueled 
combustion device are extremely Important to the performance, stability, and 
pollutant formation of the combustor. The spray characteristics should be 
measured accurately, and a model should be developed to predict the variation 
of spray characteristics, droplet trajectory, spray dispersion, and evapora- 
tion history. It Is also Important to understand the Influence of various 
factors on the measurement of spray characteristics. A comparison between 
predicted and measured spray characteristics provides useful Insight Into 
spray combustion. Such research has a significant effect on engineering 
applications and on the fundamental understanding of spray combustion. 

The authors analyzed the variation of spray characteristics along the 
axial distance downstream of a pressure-swirl atomizer. The effect of spray 
evaporation on the variation of spray characteristics Is discussed In 
reference 1 and was presented by J.S. Chin and O.Y. Zhu at this symposium. 
(Paper, entitled The Interdependence of Spray Characteristics and Evaporation 
History of Fuel Sprays In High Temperature Airflows, was unavailable for 
printing at the time of publication.) The Influence of downstream distance on 
the spray characteristics of pressure-swirl atomizers Is caused by the effects 
of spray dispersion and drop acceleration or deceleration. This condition Is 
analyzed In reference 2. These analyses have been partially substantiated by 
some experimental data. Until now there have been few experimental data of 
the variation of spray characteristics along the radial distance downstream of 
the pressure-swirl atomizer. Lee Oodge published his measurement results as 
shown In figure 1 (ref. 3). He shows that both the Sauter mean diameter SMD 
and drop-size distribution parameter N Increase with radial distance. There 
Is no analysis on this aspect available at the present time. 

We are Interested In comparing the experimental and predicted values of 
the variation of spray characteristics along the radial distance. Therefore, 
we carefully measured this variation In a well-defined experimental condition, 
developed an analytical model to predict the variation, and then compared the 
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experimental and predicted data. These comparisons promote our understanding 
of spray dispersion and the factors Influencing drop-size measurement. After 
we have predicted the simple spray system, we will be able to expand and 
modify the proposed spray model to extend Its capability to more complicated 
spray systems. 

Cao and Chin proposed a flat-fan spray model for the fuel distribution 
downstream of a plain orifice Injector under cross airflow (ref. 4). The 
model was validated by fuel distribution measured by gas analysis and drop 
size. It Is obvious that, for the validation of a given spray model, the 
spray characteristics measurement Is reliably accurate and much simpler than 
the fuel distribution measurement. It Is this author's Intention to validate 
a spray model by measuring the spray characteristics along the radial distance. 
The same spray model has also been validated by measuring the spray character- 
istics along different e angles (as shown In fig. 5). 


TEST APPARATUS, ATOMIZER, AND PARTICLE SIZER 

The test section was a rectangular chamber with transparent windows for 
optical drop-size measurement. The test apparatus Included an air system and 
a fuel system, as shown In figure 2. 

The atomizer tested (fig. 3) was a pressure-swirl atomizer taken from an 
existing aircraft gas-turbine engine. The atomizer had a pilot-fuel flow and 
a main-fuel flow. In this test we used the pilot-fuel passage. Because the 
flow rate Is relatively lower with pilot-fuel flow, the obscuration was In the 
right range. 


A Malvern drop-size analyzer model 2200 was used. The principle of the 
particle sizer Is shown In figure 4. The Instrument was set on the Rosln- 
Rammler distribution mode. The direct results from the Instrument are the 
characteristic diameter B and the drop-size distribution parameter N In 
the R-R distribution which Is expressed by 


Q = 1 



(1) 


where Q Is the fraction of the total volume contained In drops that have a 
diameter less than D. 


The Sauter mean diameter SMD can be obtained from 


SMD 



( 2 ) 


where r Is the gamma function. 


EXPERIMENTAL RESULTS 

During drop-size measurement tests, the laser light beam was first 
adjusted so that It was In the same horizontal plane as the central line of 
the atomizer and that It was perpendicular to this central line. At this 
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position (defined as the y = 0 position), the spray characteristics were 
measured. Then the Malvern particle sizer was put on some size blocks of 
known dimension. This changed the height of the laser beam so that the radial 
distance where the particle sizer measured the drop size was changed. The 
summation of the size block was defined as the radial distance y. The dimen- 
sions of the blocks chosen were 2.5, 3.0, and 10 mm. Thus, the following 
radial distance y could be obtained by different combinations of the size 
blocks: y = 0, 5.5, 10.0, 15.5, 20.0, 25.5, and 30.0 mm. The experimental 

results are shown In figures 6 to 13. From these figures It Is clear that the 
Sauter mean diameter and drop-size distribution parameter Increase with radial 
distance until the edge of the spray has been reached. The experimental 
results obtained by the present authors are quantitatively In good agreement 
with Lee Dodge's results. The obvious explanation of the measured results 
shown In figures 6 to 13 Is that with larger radial distances the possibility 
of the laser light beam meeting the large drops Is greater than It Is with 
smaller radial distances. This shows that when we measure the drop size It Is 
necessary to take a sample that Is representative of the spray formed; how- 
ever, It Is Impossible to define the sample which would fully represent the 
spray. The most reasonable sample can be obtained by measuring the spray at 
the y a 0 position. 

Comparing figure 6 with 8 and 7 with 9 shows that, at a higher airflow 
velocity, the change of spray characteristics along the radial distance Is 
weaker. A comparison of figures 8, 10, and 12 and of 9, 11, and 13 shows 
that, with a lower pressure drop across the nozzle, the change of the spray 
characteristics along the radial distance Is weaker. These results tell us 
that when the spray Is expanded less radially the change of spray characteris- 
tics Is also flatter. The reverse Is also true. 

The experimental data obtained In this way were used to check the pro- 
posed spray model . 


PHYSICAL MODEL 

The spray formed by a pressure-swirl atomizer Is shown In figure 5. The 
physical model was proposed with the following assumptions: 

(1) The atomization process Is completed as soon as the fuel leaves the 
atomizer nozzle; that Is, droplets of different sizes start their movement at 
the nozzle exit. Droplets of different sizes have the same Initial velocity. 
Because It Is suitable to use a cylindrical coordinate system to describe 
droplet motion, the tangential velocity and the radial velocity component can 
be combined and treated as one velocity component. 

(2) The airflow Is uniform. 

(3) The air Is at ambient temperature, so droplet evaporation can be 
neglected. 

(4) The Influence of different droplet velocities on the drop-size meas- 
urements made by the Malvern particle sizer Is neglected. 

The Initial drop-size distribution can be expressed by D 0 (SMD 0 ) and N 0 
with the Initial volume fraction (dQ/dD) 0 as 
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At a downstream distance x, the droplets of diameter D move from the origin 
and then distribute themselves on a circumference of radius R(D) while the 
laser light beam of diameter dj only meets some of these droplets: 


2(<P 0 ~ ‘Pi ) 


/dfl\ ‘^2 ~ *1 

\d D ) 2ir 


where 


*1 


= arcsln 


y o - I d 1 
R(D) 


1 


*2 * arCS ’ n ^ R(0) dj 


The definitions of q> 2 and <p-| are shown In figure 5. 

y 0 - 5 d j i R ' D > i y 0 * 5 d j 


then 


When 


(4) 


(5) 

( 6 ) 


<p 2 


2 


where y Q Is the radial distance and dj the diameter of the light beam. 

The summation of the liquid fraction of the droplets of various diameters 
passing through the laser light beam Is given by 



(7) 


The minimum diameter D m ^ n that the laser light beam can see Is determined 
by 


R ( D m1n> 



That Is, at the downstream distance x, If R(D) < y 0 - 1/2 dj, this droplet 
will not pass through the laser light beam. The relationship R(D) can be 
determined by the droplet motion equation. The maximum diameter that may 
exist In the spray, can be approximated by 


D max = 3.5(SMD) 


If we take the droplets within the laser light beam as a new spray (a 
sample), the liquid volume fraction of droplets of diameter D In this new 
spray (sample) Is 



Since this Is actually the drop-size distribution of the sample (new spray) 
the laser light beam Is to measure, we may obtain Q‘ as a function of D. 
Then, If we assume that the Rosln-Rammler (R-R) distribution can be used for 
data fitting the sample (new spray), we obtain new SMD* and N' for the 
sample. These values, SMD 1 and N‘, change with the distance y 0 . If we 
delete the prime, we obtain the functional relationships SMD = f(y) and 
N = f(y) . 


DROPLET TRAJECTORY EQUATION 


If we use a coordinate system which Is moving with the airflow, then the 
droplet motion equation Is 


2 

dw _ r p a 
dt D 2 


where 

3 

m droplet mass, (ir/6)p 4 D 
liquid density 

w droplet relative velocity 

Cp droplet drag coefficient, Cp = 15/^Re for Re = p a DW/v fl 
p. density of air 

a 

v- viscosity of air 

a 

Equation (9) can be rewritten as 

d ( W o) _ 45 (p a v a w o ) 1/2 n -3/2/w \ 3/2 

dt = "4 p a V“J 

where W Q Is the droplet Initial relative velocity. 


(9) 


( 10 ) 
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Integrating equation (10) yields 


S-- < Bt 


+ 1) 


-2 


(ID 


where 


/ ui \ 1/2 

n 45 ( p a v a W o) .-3/2 

B = 8 Z ° 


( 12 ) 


We know that, at the downstream distance x, droplet velocity, air veloc- 
ity, and the relative velocity of these two have the following relations: 


and 


V = V + (V - v )n- 
x a v xo a'W 

o 


V = v — 
R Ro W 

o 


(13) 


where V x and Vr are the axial and radial velocity components of the drop- 
let, and V xo and Vr 0 are the Initial values of these components. Thus, we 
obtain the droplet trajectory equations 


and 



From equation (11) we have 


(14) 


f 


W_ . t 

W dt " Bt + 1 
0 


(15) 


From equations (14) and (15) we have the radial position (droplet size 0) at 
downstream distance x: 


R(D) = V 


1 


Ro 


B + t 


(16) 
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Time t can be determined by 


XB - V 


t = 


xo 


* [(XB - V xo ) 2 ♦ 4V a Bx] 1/2 


XO' 

2V B 

a 


(17) 


where B Is determined by equation (12). 

From equations (14) and (15), by substituting x and R * = R(D * ) = 

y Q - (1/2) dj, we obtain 


min 


45 (V, M ,) 
8 p. 


1/2 


V - V » 
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(18) 


CALCULATION METHOD AND RESULTS 

In order to calculate the variation of spray characteristics along the 
radial distance the following parameters must be known: 

(1) Axial distance x 

(2) Pressure drop across the nozzle APj, 

(3) Fuel properties and physical properties of air 

(4) Atomizer spray cone angle 

(5) Initial spray characteristics (SMD) 0 and N 0 

The calculation method used In this paper may be summarized as follows: 

(1) Use the measured SMD and N values at y = 0 and the downstream 
distance x as the approximation of (SMD) 0 and N 0 . We know that there Is 
some difference between (SMD) 0 and SMDy_o, N 0 , and Ny_g (ref. 2). So far 
these are the most reasonable approximate values that can be used. 

(2) Calculate the Initial droplet velocity Vj, = y^Apg/pj.)^ where y 
Is the flow coefficient for the pressure drop. For the atomizer used In the 
present research, the spray angle Is close to 90°; therefore, 


( 19 ) 


(3) Determine the D m ^ n value for a different y 0 value by 
equation (18). We assume D^ax * 3 . 5( SMD) . 

(4) Use the Simpson method to solve equation (8), with the R(D) value 
calculated from equations (16) and (17). 
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(5) Determine the Integration 



From this, we are able to obtain a set of the data of Q' and D. 
(6) Assume that Q' can be fitted by R-R distribution: 



By using the least-square method, we can obtain D and N ' ; then, SMD' and 
N 1 can be calculated. These SMD 1 and N 1 values are plotted against y 
In figures 6 to 13. 

The results of the calculations for different airflow velocities and for 
the nozzle pressure drop are shown In figures 6 to 13. These are compared 
with the corresponding experimental data. It Is obvious that the predictions 
are In excellent agreement with the experimental results. 


CONCLUDING REMARKS 

In the present study, droplet evaporation was neglected (discussed by 
J.S. Chin and J.Y. Zhu In the unpublished paper cited previously), and It Is 
not difficult to take Into consideration the effect of droplet evaporation. 

After this model has been validated, It Is relatively simple to predict 
the fuel distribution downstream of a pressure-swirl atomizer (ref. 6). 

One Interesting point that was noticed from the experiments and the anal 
ysls was that the larger droplets are always at the spray edge. Because the 
Ignition spark plug Is also always positioned close to the spray edge, the 
spark might meet the larger droplets, which would harm the Ignition. 

We measured the variation of spray characteristics along the radial dis- 
tance was by laser light-scattering technology. And we found that the predlc 
tlon based on the proposed spray model was In excellent agreement with the 
experimental data. 
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Figure 2. - Experimental system. Plenum chamber, 1; air 
temperature, 2; air flowmeter, 3; test section, 4; exhaust 
pipe, 5; fuel tank, 6; filter, 7; pump, 8; return valve, 9; 
valve, 10; pilot-fuel flow valve, 11; pressure gauge, 12; 
atomizer, 13; main-fuel flow valve, 14; pressure gauge, 15. 





Figure 3. - Atomizer. Pilot-fuel 
flow, 1; main-fuel flow, 2. 
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Figure 4. - Laser drop particle sizer. Laser, 1; beam expander, 2; test 
section, 3; fuel nozzle, 4; lens, 5; photo detector, 6; solid valve, 7; 
amplifier, 8; A/0 converter, 9; data transfer, 10; Interface, 11; 
computer, 12; printer, 13; micro tune, 14. 













Figure 6. - Variation of Sauter mean diameter (SMD) along radial 
distance (y) for x = 70 mm, V a = 9 m/sec, and AP = 12 kg/cm?. 



Figure 7. - Variation of drop-size distribution parameter (N) 
along radial distance (y) for x = 70 mm, V a = 9 m/sec, and 
AP = 12 kg/cm2. 
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Figure 10. - Variation of Sauter mean diameter (SMD) along radial 
distance (y) for x = 70 mm, V a = 6 m/sec, and AP = 8 kg/cm2. 
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Figure 11. - Variation of drop-size distribution parameter (N) 
along radial distance (y) for x = 70 mm, V a = 6 m/sec, and 
AP = 8 kg/cm2. 
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Figure 12. - Variation of Sauter mean diameter (SMD) along radial 
distance (y) for x = 70 mm, V a = 6 m/sec, and AP = 4 kg/cm 2 . 
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Figure 13. - Variation of drop-size distribution parameter (N) 
along radial distance (y) for x = 70 mm, V a = 6 m/sec, and 
AP = 4 kg/cm 2 . 
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An experimental program to characterize the spray from candidate nozzles 
for Icing-cloud simulation Is discussed. One candidate nozzle, which Is cur- 
rently used for Icing research, has been characterized for flow and drop size. 
The median-volume diameter (MVD) from this air-assist nozzle Is compared with 
correlations In the literature. The new experimental spray facility Is dis- 
cussed, and the drop-size Instruments are discussed In detail. Since there Is 
no absolute standard for drop-size measurements and there are other limita- 
tions, such as drop-size range and velocity range, several Instruments are used 
and the results are compared. The drop-size Instruments used are the Malvern 
2600, the NASA Lewis Scanning Radiometer, the Phase/Doppler Spray Analyzer 
(P/DSA), a Forward-Scattering Spectrometer Probe (FSSP), an Optical Array Probe 
(OAP) , and a limited amount of photography. The known capabilities and limi- 
tations of the drop-size Instruments are compared. 

A two-phase model was developed at Pennsylvania State University. The 
model uses the k-e model of turbulence In the contlnous phase. Three methods 
for treating the discrete phase are used: (1) a locally homogenous flow (LHF) 

model, (2) a deterministic separated flow (OSF) model, and (3) a stochastic 
separated flow (SSF) model. In the LHF model both phases have the same veloc- 
ity and temperature at each point. The DSF model provides Interphase transport 
but Ignores the effects of turbulent fluctuations. In the SSF model the drops 
Interact with turbulent eddies whose properties are determined by the k-e 
turbulence model. The two-phase flow model has been extended to Include the 
effects of evaporation and combustion. Model development has been accompanied 
by an experimental program to provide data for verification. Presently, the 
model has been further extended at NASA Lewis to Include the effects of 
swirling flow, and a particle-laden swirling flow experiment Is currently 
being performed to provide verification data. 


INTRODUCTION 

Sprays and two-phase flow are of Interest In many fields - for example, 
gas turbine combustion, agricultural crop spraying, paint spraying, fire pro- 
tection, and drying. This paper describes the use of sprays to simulate Icing 
clouds for the Altitude Wind Tunnel Icing Research Program. Icing-cloud simu- 
lation and gas turbine combustion research have common Interests In the method 
of Injection, the characterization of drop size and of droplet velocities, and 
the development of analytical models for droplet trajectories and droplet 
vaporization. This paper discusses Injector or spray-nozzle drop-size charac- 
terization, an experimental test facility for Injector characterization, drop- 
size measurement techniques, and analytical models for trajectory calculations. 


INJECTOR CHARACTERIZATION 

The Injector, or spray nozzle, used currently for Icing-cloud simulation 
at the NASA Lewis Research Center Is an air-assist Injector (fig. 1). It Is 
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very similar to the Injector used by Nuklyama and Tanasawa (ref. 1). The water 
flows through a center tube (0.635-mm dlam) and has a concentric annular pas- 
sage (9.35-mm dlam) for air-assist flow. The water and the air exit through a 
3.175-mm orifice. The advantage of an air-assist Injector for Icing research 
and other applications Is that, within limits, there Is nearly Independent 
control of water flowrate and drop size. The water flow Is controlled by the 
difference between the air pressure and water pressure, and the drop size Is 
controlled by the air pressure only. 

Drop-size measurements using this Injector were made for us at Arnold 
Engineering Development Center by J.L. Hunt. The drop-size data shown In 
figure 2 are plotted as median-volume diameter In micrometers. The median- 
volume diameter MVD Is defined as the drop size where 50 percent of the total 
droplet volume Is contained In droplets less than that diameter, and 50 percent 
Is contained In droplets greater than that diameter. The MVD Is plotted as a 
function of the difference of water and air-assist pressure with the air- 
assist pressure as a parameter. 

This method of correlating the data Is very useful for setting the oper- 
ating conditions In the wind tunnel. However, the drop-size data could be 
presented In a more fundamental manner as a function of the water to air mass 
ratio. Other factors that appear In correlations In the literature, such as 
viscosity, surface tension, and relative velocity between the air and water, 
are constant or nearly constant for these data. The relative velocity Is 
nearly constant, even though the flowrates are changing. The reason for this 
Is, first, that the air Is being Injected at high pressure Into a low-pressure 
ambient environment and thus has sonic, or nearly sonic, velocity. Second, the 
water pressure Is relatively low so that the water has a relatively small 
velocity compared to the sonic velocity of the air. Thus, the difference In 
the velocities or the relative velocity Is constant. 

These data were compared to values from various correlations. These com- 
parisons are shown In figure 3 for one air-assist pressure. The Injectors used 
In the references were similar to the NASA Lewis Research Center Injector. 
However, the drop-size measurement techniques were different. The Nuklyama- 
Tanasawa (ref. 1) drop-size correlation was made from measurements of water 
droplets on photographs. The Lorenzetto-Lef ebvre correlation (ref. 2) was made 
from light-scattering measurements of water droplets. Wlgg (ref. 3) used 
measurements from molten wax sprays with swirl. As you can see In the figure, 
the results of the correlations have a wide range of values, some higher and 
some lower than the data taken at Arnold. Preliminary data taken at NASA 
Lewis seem to be closest to the data of Nuklyama-Tanasawa. Preliminary 
data taken with the NASA Lewis Injector at Carnegle-Mellon University are lower 
than the Arnold data. There needs to be much more work on drop-size measure- 
ments and correlations before drop sizes can be accurately predicted. 


TEST FACILITY 

The objectives of our experimental program are to characterize a variety 
of candidate Injectors for drop size and dispersion, to compare different 
drop-size measurement techniques, and to take data that can be used for tra- 
jectory models. Trajectory model data consist of Initial condition data taken 
upstream close to the Injector and verification data taken far downstream. 
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The test stand Is shown In figure 4. Duct air enters through a 40-cm- 
dlameter pipe. The maximum air velocity through the pipe Is 88 m/sec. We will 
be able to control the relative humidity from 0 to 100 percent by adding steam 
far upstream. Controlling the relative humidity could be necessary If evap- 
oration Is Important. 

The duct air temperature will be controlled by Injecting liquid nitrogen 
upstream. The primary use of the temperature control will be to provide a 
constant Inlet temperature of 292 K. For some Icing tests the Inlet temper- 
ature Is held at 273 K and even lower temperatures for shorter times. 

The Injector Is Inserted In a section with a window so that the nozzle 
operation can be monitored with a television camera. Water Is the only liquid 
capable of being tested In this facility. The water Is demineralized and 
heated to 380 K. For Icing applications, the water and the air-assist Injector 
are heated so that Ice crystals do not form during the expansion process. The 
water and air-assist flow temperature can be controlled between the ambient 
temperature and 373 K. 

The Injector Is shown In a section 2.5 m upstream of the Instrument sec- 
tion. The Injection section can be moved to any location closer to the 
Instrument section. Upstream of the Injector there Is a traversing pitot tube 
to measure the Inlet velocity. 

Downstream, at the Instrument section, drop-size measurements can be made. 
The Instrument section has a square cross section with windows In the sides 
that are at 7° angles to each other to prevent reflections. Other drop-size 
Instruments that are used are described In the next section. The Instruments 
sit on a three-axis traversing table so that radial and spatial variations can 
be measured. A liquid-water measurement probe traverses vertically to measure 
the variation In liquid-water concentration. 

Another set of droplet-sizing Instruments Is mounted In the duct down- 
stream of the Instrument section. The principle of operation of these Instru- 
ments Is discussed In the next section. These Instruments are Intrusive and 
mounted In a fixed position. They are the same Instruments that are mounted 
on airplanes to measure drop sizes In clouds. Tests with these Instruments 
will provide a good comparison of the data from the nozzles with that In the 
clouds. This Is an atmospheric pressure rig, so the air and water are vented 
to the atmosphere. 


DROP-SIZE INSTRUMENTS 

The most Important factor In making drop-size measurements Is the drop- 
size measuring Instrument. As you have seen from the figure comparing the 
drop-size data correlations, there must be large differences In the data on 
which they are based. One of the major problems Is that there Is no standard 
to which to compare the drop-size measurements, especially In a spray. Since 
there are other limitations beside the lack of a standard such as drop-size 
range and velocity range, In our program we use five different drop-size 
Instruments and compare the results. The Instruments are a Malvern 2600, a 
NASA Lewis Scattered-Light Scanner, a Phase/Doppler Spray Analyzer, a Forward- 
Scattering Spectrometer Probe, and an Optical Array Probe. And a limited 
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amount of photography Is used for qualitative analysis. A qualitative analysis 
and a discussion of each of these Instruments follows. 

The Malvern 2600 Is a commonly used drop-size Instrument (fig. 5). Lewis' 
HSD version of this Instrument, which provides a faster snapshot of the data. 

Is described In reference 4. (The Inserts In figure 5 describe the principle 
of operation.) This Instrument measures the angle that the light Is diffracted 
by the drops. The bigger drops diffract the light at a small angle and the 
smaller drops diffract the light at a larger angle. The detector Is a series 
of concentric rings that determine the angle of the scattered light. The drop 
size Is calculated from the scattered light angle. The Insert also shows that, 
no matter where the droplet Is In the line of sight, the light that Is scat- 
tered at the same angle always falls on the same detector ring. 

Although there Is no standard for determining drop sizes, we use a 
Hlrleman reticle for calibration purposes. The Hlrleman reticle (fig. 6) con- 
sists of a glass substrate on which chrome disks are photographically depos- 
ited. The chrome disks are sized to represent a Rosln-Rammler distribution 
with an x-bar of 50 and an n value of 3. One problem we have found using 
the reticle with the Malvern 2600 Is that there Is a correction for anomalous 
diffraction. The Malvern calculations are based primarily on Fraunfaufer dif- 
fraction theory. However, for small drop sizes (below a diameter of 5 ym) 
three-dimensional effects are considered; these effects are lumped under the 
term anomalous diffraction. Thus, the problem Is that the reticle has two- 
dimensional disks and not three-dimensional spheres. We are currently studying 
the feasibility of modifying the Malvern software to Ignore the anomalous dif- 
fraction calculation when using the reticle. If this does not prove to be 
practical, we will still use the reticle for relative comparisons. 

We also plan to add an additional lens system to the Malvern to be able 
to measure an MVD of 5 ym across the 30-cm Instrument section. Verification 
of the correctness of the additional lens system will be done by taking meas- 
urements of the reticle with and without the additional lens system. 

The second Instrument, the Scattered-Light Scanner, Is also based on the 
diffraction principle. This Is a revised version of the NASA Lewis Scanning 
Radiometer (described In ref. 5), which was developed by Don Buchele and used 
In many research projects by Robert Ingebo. Figures 7 to 9 explain the oper- 
ation of the Instrument. 

Figure 7 shows the reflection, refraction, and diffraction of a laser 
light beam passing through a drop. Since calculations are based on the dif- 
fraction of light through the drop, the wide-angle diffraction Is Ignored. 

There Is narrow-angle diffracted light from the envelope (or aperture) of the 
laser beam that must also be Ignored. 

Figure 8(a) Is a plot of light Intensity as a function of the light-scat- 
ter angle. Figure 8(b) shows the product of the light Intensity and the light- 
scatter angle as a function of the light-scatter angle. It shows that the 
measurement Interval extends from where the Intensity of the near envelope 
diffraction Is negligible to a value of the light Intensity which Is still high 
enough to be meaningful. A scanner Is used to measure the Intensity of the 
scattered light as a function of the light scatter angle. The scanner Is a 
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disk with a radial slot that Is rotated with Its centerline eccentric to the 
laser beam centerline. 

To calculate the mean drop size, Don Buchele has calculated the light 
Intensity as a function of a beam spread parameter (fig. 9). From the previous 
figure we can find the maximum light Intensity and the light Intensity at a 
specified angle. Using these values and the beam spread parameter and knowing 
the wave length and the specified angle, we can determine the mean drop size. 
Calculations were made with five different distribution correlations, and they 
are shown plotted for three dispersion levels. The distribution correlations 
used were the Nuklyama-Tanasawa, the upper-limit distribution functions, the 
Rosln-Rammler, the Welbel, and the Rinks distribution function. The beam- 
spread parameter Is nearly the same for all the distributions and dispersions 
at a relative Intensity of 0.75. If this relative Intensity Is used, then the 
mean drop-size calculation Is Independent of distribution and dispersion. The 
dispersion can be calculated from the slope of the curve. 

A schematic of the Scattered Light Scanner (fig. 10) shows the various 
parts and dimensions of the Instrument. The primary feature Is the eccentric 
light scanner. 

The third drop-size Instrument to be discussed was developed under a NASA 
combustion fundamentals contract to William Bachello at Aerometrlcs Corpora- 
tion. The method Is the phase/Doppler detection technique and the Instrument 
Is called the Phase/Doppler Spray Analyzer (ref. 6). This Is a scheme where a 
set of laser beams Is crossed to form a fringe pattern (fig. 11). As a drop- 
let passes through the fringes the light Is scattered and measured 30° off axis 
with a photodetector. From this measurement, simultaneous values of drop 
diameter and velocity can be calculated. The measurement and calculations are 
for single droplets; averages and medians are calculated from the single drop- 
let data. 

Figure 12 shows In more detail the principle of the phase/Doppler system. 
The scattered light Is collected by three detectors with slightly different 
spatial locations. Because of the spatial difference, the light path to the 
three detectors Is different, and thus there Is a phase difference In the sig- 
nals. From this phase difference the drop size can be calculated. Again, the 
phase difference occurs because of the path difference Inside the drop, so that 
the phase angle Is a function of the diameter. 

Figure 12(b) Is a plot of phase difference as a function of diameter. 

From the phase difference between detectors 1 and 2, an approximate drop-size 
estimate can be made. From this approximate drop size and the phase difference 
between detectors 1 and 3, a more accurate determination of diameter can be 
made. 


The calculations are based on reflection and refraction theory. Diffrac- 
tion can be Ignored because the measurements are taken 30° off axis, where 
diffraction Is negligible. 

The Forward-Scattering Spectrometer Probe (FSSP) Is one of the two Intru- 
sive Instruments that we use that Is located downstream of the Instrument sec- 
tion. This Instrument Is the same type that Is mounted on airplanes to measure 
raindrop diameters (see ref. 7). It provides a comparison of the other drop- 
size Instrument measurements with the measurements made In clouds. 
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The FSSP Is a single particle counter. A schematic of Its operation Is 
shown In figure 13. A laser beam. Is optically directed across two arms of the 
probe, In a direction perpendicular to that of the particle-laden flow. Scat- 
tered light from a particle Is optically directed and focused onto a photodiode 
to measure the Intensity of the scattered light. Any unscattered light Is 
removed by the dump spot. The scattered light Is also directed onto a second 
photodetector that Is used to determine whether the droplet Is In focus (In 
focus along the axial length of the laser beam). This second photodetector 
also has a dump spot. The light from a droplet that is In focus Is removed by 
the dump spot. A droplet that Is not In focus has a broader signal that Is not 
removed by the dump spot, and the droplet measurement Is rejected. 

Drop size as a function of the measured Intensity Is shown In figure 14. 

In the figure the light Intensity Is shown In the near forward-scattering 
angles of 4° to 14°. The curve Is shown for a Gaussian T-type (transverse 
electric) beam. This Is a multivalued function; thus, for a given Intensity, 
the drop size may not be uniquely determined. The actual Instrument uses a 
high-mode laser beam which has a smooth curve and Is single valued. 

The last drop-size Instrument to be discussed Is the Optical Array Probe 
(OAP) , which can also be flown on an airplane to measure raindrop sizes. This 
Instrument Is made by the same company as the FSSP, but It Is made to cover a 
larger drop-size range. The operating principle Is shown In figure 15. It Is 
a single-particle Imaging system. The drops pass through a laser beam, and 
then the light Is transmitted through the Instrument optics to an array of 
photodetectors. The drop size Is determined by counting the detectors In the 
shadow of the drop. The last detector at each end Is kept open as a check to 
determine If the whole drop Is seen. To determine If the drop Is In focus, we 
use an Intensity gradient criteria. An Infocus drop has a sharp gradient from 
lit to unlit photodetectors, while an out-of-focus drop has a gradual gradient 
from lit to unlit. Figure 16 Is a simple schematic of the OAP. 

Table I compares the FSSP, OAP, Malvern 2600, Scattered-Light Scanner, and 
Phase/Doppler Particle Analyzer (P/DSA) drop-size Instruments. The size range 
refers to the droplet diameters that can be measured with any particular setup. 
The number of bins refers to the number of Increments Into which the size range 
Is divided. The overall size range then refers to the smallest to the largest 
diameter that can be measured, but It may require a change In setup and multi- 
ple measurements to cover the range. The numbers In the ranges refer to Indi- 
vidual droplets for all the Instruments except for the Scattered Light Scanner, 
whose range Is given In MVO. 

The velocity limits are all acceptable for our work except for those of 
the P/DSA, which has a slight limitation for our purposes. Currently manufac- 
tured P/DSA Instruments have a upper limit of 75 m/sec. The comments section 
lists features that we consider particularly beneficial. With the FSSP and 
OAP, we consider the ability to use the same Instrument to obtain data from 
clouds and the test facility to be Important. The Malvern and the Scattered- 
Light Scanner have the ability to measure MVD directly and thus are good for 
screening a number of Injectors quickly. The P/DSA has the capability of 
making point measurements of size and velocity which makes the data particu- 
larly useful for modeling work. 


68 


ANALYTICAL TRAJECTORY MODELS 


The second part of this report describes spray model development. The 
model work to he presented was done at Pennsylvania State University by Pro- 
fessor Faeth and his graduate students and associates (refs. 8 to 10). 

Figure 17 shows the Important Interactions of the drop with Its environment. 
These Interactions are: the droplet drag that results In a source term In the 

continuous flow equation, random deflections of the drops by turbulent eddies, 
and turbulence modulation of the continuous phase by the extraction of turbu- 
lent energy by the droplet drag. Actually, three models of Increasing com- 
plexity were developed. They are a locally homogenous flow model, a determin- 
istic separated-f low model, and a stochastic separated-flow model. 

The locally homogenous flow model Is the simplest model (fig. 18). It 
uses the k-e model to calculate the continuous-phase flow characteristics. 

It then assumes that the particles have the same velocity as the continuous 
phase and treats the k-e equations as a varylng-denslty flow. This model Is 
correct for small drops that follow the continuous phase flow. An advantage 
of the model Is that the calculation Is no more complex for a two-phase flow 
than for a single-phase flow. It Is also useful when the Initial conditions 
of the spray, such as drop size or spray angle, are not available since they 
are not needed In the calculation. 

The second model Is the deterministic separated flow (DSF) model 
(fig. 19). It also uses the k-e model to calculate the gas or continuous- 
phase flowfleld. The particle or droplet distribution Is divided Into a number 
of classes with a representative number of droplets In each class. Then, for 
each drop, the trajectory Is calculated. The trajectory calculation Is based 
on the mean continuous-phase flow properties. Dispersion due to turbulence Is 
Ignored. This Is the most widely used approximation In current models of com- 
busting sprays. Other assumptions are that the drops are spherical and that 
the spray Is dilute with no collisions. 

One feature of this model Is the partlcle-ln-cell approach. As a particle 
goes through a cell, the drag momentum Is coupled with the continuous-phase 
momentum. The net change In momentum Is calculated for each drop class as It 
passes through a computational cell as shown In figure 20. Then the gas or 
continuous-phase flow equation Is calculated with the change In momentum as a 
source term. 

The last, and most complex, model Is the stochastic separated flow (SSF) 
model (fig. 21). In this model, the continuous-phase turbulence does Interact 
with the drop or particle. Again the gas-phase calculation Is made using the 
k-e model. To make the Lagraglan drop-trajectory calculation we assume the 
drop Interacts with an eddy. The properties of the eddy have a distribution 
based on a probability density function with a standard deviation dependence 
on k. From the distribution, the properties of the eddy are found by using a 
Monte Carlo technique. The Interaction Is shown In figure 22. The drop 
Interacts with the eddy either for the lifetime of the eddy or the time It 
takes the particle to traverse the eddy. The size of the eddy Is based on the 
dissipation length. The number of calculations are much greater for this 
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method than for the deterministic separated flow (DSF) model. While for the 
DSF model approximately 600 drops are sufficient for an accurate calculation, 
for the stochastic separated flow model about 6000 drops are needed for an 
accurate calculation. 

An experimental study was also done to complement the analytical modeling 
effort. From a large amount of data, four plots were chosen as representl ve. 
Figure 23 Is a plot of the mean axial velocities of the drops as a function of 
the drop diameter at different axial locations. An air-assist Injector with 
an exit-orifice diameter of 1.194 mm was used. The Initial conditions were 
taken at an x/d of 50, where x Is the axial distance and d Is the Injec- 
tor-orifice diameter. This Is a distance of about 55 mm. This was sufficient 
distance for the spray to be well formed for most conditions. However, there 
were some conditions where ligaments were observed. 

The drop-size and velocity measurements were made using double-flash 
photography. Two other methods were used to measure drop size. One was slide 
Impaction In which a slide coated with magnesium oxide was used. The second 
alternate method was to measure drop size with a Malvern. The Malvern meas- 
urement was used mainly to monitor the spray to determine whether the spray 
varied from day to day. 

The data were compared with the calculations from the stochastic separated 
flow model and the deterministic separated flow model. The calculations agree 
well, with the data from the stochastic separated flow model agreeing the best. 
For all the data plots, the stochastic separated flow model agreed well with 
the data, while the deterministic separated flow model and the locally homoge- 
nous model agreed well for some conditions and not so well for other condi- 
tions . 

Figure 24 shows the radial variation of mean liquid flux 6. For this 
case, the SMD at an x/d of 50 was 30. Again the calculations with the 
stochastic separated flow model agree well with the data for all conditions, 
while the other two models do not agree well for all conditions. It Is Inter- 
esting to note that the locally homogenous model underpredicts the data. For 
all the previous work In which the locally homogenous model was calibrated with 
sol Id-particle-laden flow, the model overpredicted the development of the 
spray. 

Figure 25 shows the variation of the Sauter mean diameter with the dis- 
tance along the axis. Figure 25(a) Is for a nonevaporating oil spray and 
figure 25(b) Is for an evaporating Freon 11 spray. Case 1 had an Initial SMD 
of 30 at an x/d of 50 and case 2 had an Initial SMD of 87. For the nonevap- 
orating case, the mean diameter Increases after an x/d of 100 because of the 
greater dispersion of the small drops. The stochastic separated flow model 
does better at predicting the data because It considers the dispersion due to 
turbul 9HC 9 . 

The variation of SMD along the axis for an evaporating spray Is shown In 
figure 25(b). For an evaporating spray, there are competing effects: evapo- 

ration reduces the mean diameter, and greater dispersion of the small drops 
Increases the mean diameter along the axis. The stochastic separated flow 
model and deterministic separated flow model do a good job of predicting the 
trends even If they do not exactly duplicate the data. 
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The model development Is continuing at the NASA Lewis Research Center as 
part of a Ph.D. thesis (by Daniel Bulzan). The ability to calculate flow 
fields with swirl has been added to the model, and an experimental calibration 
and verification Is underway. The experiment will have a weakly swirling flow 
field with a swirl number of less than 0.4. Glass beads with known sizes will 
be used to simplify the analysis. Particle number flux will be measured by a 
30° off-axis Mle-scatterlng technique. Laser Doppler veloclmetry (LDV) will 
be used to measure the particle and gas-phase velocities. The parameters to 
be studied are swirl number, particle loading ratio, and particle size. The 
Initial conditions will be well characterized to provide Input conditions to 
the computer model. 

Figure 26 Is a drawing of the experimental configuration. It has a two- 
component laser veloclmeter system that will be used to measure the particle 
and the gas-phase velocities. Small aluminum oxide particles, approximately 

0.5 ym diameter, will be added to make the gas-phase velocity measurement. The 
laser system will remain fixed and the Injection system will be traversed to 
determine spatial variations. The Injector has four tangential slots to create 
swirl . 


CONCLUDING REMARKS 

There has been much work In this field of two-phase flow In the past, 
and, with the new Instruments that are more accurate, faster, and easier to 
use, there promises to be great advances In the experimental data base. Ana- 
lytical models will complement the experimental work; there promises to be more 
accurate and faster analytical models to calculate the two-phase flow fields 
that exist In Icing clouds, gas turbine combustors, and other applications. 


REFERENCES 

1. Nuklyama, Shlro; and Tanasawa, Yasushl (E. Hope, trans.): Experiments on 

the Atomization of Liquids In an Air Stream. Rept. No. 4, The Influences 
of the Characteristics of the Liquid on the Atomized Droplets. Defense 
Res. Board, Dept. Nat. Defense, Ottawa (Canada), Mar. 18, 1950. (Trans, 
from Trans. Soc. Mech. Eng. (Japan), vol. 5, no. 18, Feb. 1939, pp. 68-75.) 

2. Lorenzetto, G.E.; and Lefebvre, A.H.: Measurements of Drop Size on a 

Plain-Jet Alrblast Atomizer. AIAA Paper 76-679, July 1976. 

3. Wlgg, L.O.: Drop-size Prediction for Twin-fluid Atomizers. J. Inst. 

Fuel, vol. 37, no. 286, Nov. 1964, pp. 500-505. 

4. Watson, D.J.: Laser Diffraction Measurements In Transient Spray 

Conditions. ICLASS-85 . 3rd International Conference on Liquid Atomization 
and Spray Systems, Vol. 1, P. Elsenklam and A. Yule, eds.. Institute of 
Energy, London, 1985, pp. VC/4/1 - VC/4/15. 

5. Buchele, D.R.: Particle Sizing by Measurement of Forward-Scattered Light 

at Two Angles. NASA TP-2156, 1983. 


71 


6. Bachalo, W.D.; and Houser, M.J.: Phase/Doppler Spray Analyzer for 

Simultaneous Measurements of Drop Size and Velocity Distributions. Opt. 
Eng. . vol . 23, no. 5, Sept. -Oct. 1984, pp. 583-590. 

7. Knollenberg, R.6.: The Use of Low Power Lasers In Particle Size 

Spectrometry. Practical Applications of Low Power Lasers, SPIE vol. 92, 
D.D. Eden and J.S. Chlvlan, eds., SPIE, 1976, pp. 137-152. 

8. Shuen, J-S., et al . : Structure of Particle-Laden Jets: Measurements and 

Predictions. AIAA J., vol. 23, no. 3, Mar. 1985, pp. 396-404. 

9. Solomon, A.S.P., et al.: Measurements and Predictions of the Structure of 

Evaporating Sprays. J. Heat Trans., vol. 107, no. 3, Aug. 1985, 

pp. 679-686. 

10. Shuen, J-S.; Solomon, A.S.P.; and Faeth, G.M.: Drop-Turbulence 

Interactions In a Diffusion Flame. AIAA J., vol. 24, no. 1, Jan. 1986, 

pp. 101-108. 


TABLE I. - DROPLET-SIZE INSTRUMENT COMPARISON 



Forward-Scattering 
Spectrometer Probe, 
FSSP 

Optical Array 
Probe, 

0AP 

Malvern 2600 

Scattered-Light 

Scanner 

Phase/Doppler 
Particle Analyzer, 
P/DPA 

Size range 

1 to 16 

2 to 32 
2 to 47 
5 to 95 

20 to 300 

1.3 to 128 
1.9 to 188 
5.7 to 564 

6 to 120 MVD 

^max^^min “ 35 

Number of bins 

15 

15 

15 

NA 
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Overal 1 size range 

1 to 95 

20 to 300 

1.3 to 564 

6 to 250 MVD 

1 to 1000 

Velocity 1 imi ts 

5 to 100 m/sec 

5 to 100 m/sec 

None 

None 

3.2 x fringe 
spacing 41 m/sec 

Measurement 

theory 

Scattering 
ampl itude 

Particle 
shadow size 

Fraunhofer 

diffraction 

Fraunhofer 

diffraction 

Geometric optics 

Type of 
measurement 

Single- 

particle 

counter 

Single- 

particle 

counter 

Ensemble 

average 

Ensembl e 
average 

Single-particl e 
counter 

Probe volume 

0.218-mm (diam) 
x 2.61 mm 

4.8 mm (diam) 
x 6 cm 

Line of sight 
x 9 mm diam 

Line of sight 
x 50 mm or 
18 mm 

0.160 mm (diam) 
x 0.050 mm 

Comments 

FI iaht-instrument 

icing research 

Measures volume 
distribution 
(best for MVO) 

Size-velocity 

correlation 
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Figure 2. - Drop sizes from IRT nozzle. 
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CD-K- 1 J 6 S Figure 3. - Comparison of IRT nozzle drop-size data with 

correlations. 
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WATER — 1 1 — AIR ASSIST 
MAXIMUM PRESSURE, 1.5 MPa 
MAXIMUM TEMPERATURE, 393 K 

Figure 4. - Experimental spray facility. 



PRINCIPLE OF DIFFRACTION SIZING MEASURING MOVING PARTICLES 



SMALL PARTICLES DIFFRACT LIGHT WHEN A FOURIER TRANSFORM 
THROUGH LARGE ANGLES AND LENS IS USED, THE LIGHT DIF- 

LARGE PARTICLES THROUGH SMALL FRACTED BY PARTICLES OF A 
ANGLES GIVEN SIZE AT ANY POSITION 

IN SPACE ALWAYS FALLS ON 
THE SAME PART OF THE DE- 
TECTOR 

Figure 5.- Malvern spray analyzer. 



SIZING 

ARRAY 


Figure 6. - Photograph of sample area and quality control array of Hirleman 
calibration reticle (RR-50-3.0-0.08-102-CF-#115). Magnification lOx. 
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WIDE-ANGLE REFRACTION 
THROUGH DROP, CAN BE 
NEGLECTED 


MEDIUM-ANGLE DIFFRACTION 
AROUND DROP, 0 - x/7rD 
MEASURED 


SMALL-ANGLE DIFFRACTION 
NEAR ENVELOPE OF LASER BEAM, 
PREVENTS MEASUREMENT AT 
SMALL ANGLES 

Figure 7. - Three principal components of scattered light at test section. 
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(b) Product of light intensity and light-scatter angle 
as function of light-scatter angle. 

Figure 8. - Intensity of forward-scattered light. 
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DISPERSION - 2(D 9 - D. ^/(D 9 + D j) 



BEAM-SPREAD PARAMETER, P- JL D u 5 0 

X V - J 


Figure 9. - Envelope of scattered light for five volume-distribution 
functions and three dispersions. 
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Figure 10. - Lewis Scanning Radiometer spray analyzer. 
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Figure 16. - Optical Array P 
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2 (OAP) operating principle. 



CONDENSING 
LENS; F- 320 mm 


MIRROR 



dW 

VU MIRROR 


OBJECTIVE 


LENS; F- 60 mm 


(OAP) optical path. 




Figure 17. - Phase interactions for particle-laden flow. 



Figure 1& - Locally homogeneous flow (LHF) model. 



ASSUMPTION: FINITE INTERPHASE TRANSPORT 
RATES BUT PARTICLES ONLY 
RESPOND TO MEAN MOTION 



NOTES: I - PARTICLE DISPERSION IGNORED - ONLY 
VALID FOR "LARGE" PARTICLES 

II - MOST WIDELY USED APPROXIMATION 
IN CURRENT MODELS OF COMBUSTING 
SPRAYS AND PULVERIZED COAL 

III - FLOW FIELD FOUND FROM EULER IAN 

CALCULATION WITH DISTRIBUTED SOURCE 
TERMS DUE TO PARTICLES. LAGRANGIAN 
CALCULATION DETERMINES PARTICLE 
TRAJECTORIES 

Figure 19. - Deterministic separated flow (DSF) model. 
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ASSUMPTION: FINITE, INTERPHASE TRANSPORT RATES WITH PARTICLES 
INTERACTING WITH INDIVIDUAL EDDIES WHOSE PROPERTIES 
ARE FOUND BY RANDOM SAMPLING OF LOCAL TURBULENCE 
PROPERTIES 



+ U): 


(Up + u* p )j * (u + u')j 


NOTES : I - INVOLVES ADDITIONAL ASSUMPTIONS CONCERNING THE 
PROPERTIES, SIZE. AND LIFETIME OF EDDIES 

II - COMPUTATIONS SIMILAR TO DSF MODELS. BUT MONTE CARLO 
TECHNIQUES USED TO FIND STATISTICALLY SIGNIFICANT 
NUMBER OF PARTICLE TRAJECTORIES - CAUSING ADDED 
COMPUTATIONS 

Figure 21. - Stochastic separated flow (SSF) model. 
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Figure 22. - Particle-eddy interaction. 
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Sauter mean diameter, SMD, pm 




(b) Evaporating sprays. 

Figure 25. - Sauter mean diameter (SMD) along axis. 
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MOVABLE SCREENED ENCLOSURE 


INJECTION TUBE WITH 
AXIAL TRAVERSE 


LDA RECEIVING OPTICS 

MINICOMPUTER 


TERMINAL 



LDA SENDING OPTICS 
2-W ARGON-ION LASER 


✓ ncpii mcrnDC — > / 
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ENCLOSURE TRACK / / DADTirtc INTEGRATING DIGITAL 
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Figure 26. - Particle-laden swirling-flow experiment. 
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COMBUSTION RESEARCH ACTIVITIES AT THE GAS TURBINE RESEARCH INSTITUTE 

Shao Zhongpu 

Gas Turbine Research Institute 
Jlangyou County of Sichuan, People's Republic of China 


The Gas Turbine Research Institute (GTRI) Is responsible mainly for basic 
research In aeronautical propulsion. The combustion research activities pri- 
marily cover the following areas: annular diffuser for the turbofan augmentor, 

combustor Ignition performance, combustor airflow distribution, fuel Injectors, 
vaporizer fuel Injector, and airblast atomizer. 


RESEARCH ON THE ANNULAR DIFFUSER FOR THE TURBOFAN AUGMENTOR 

It Is well known that the Inlet velocity profile has a very strong Influ- 
ence on diffuser performance. The Inlet velocity and temperature profiles are 
not uniform: both change with flight condition. The diffuser enlarges the 

Inlet velocity distortion; that Is, the nonuniformity of the velocity Increases 
with the diffusion process occurring In the diffuser. The nonuniform Inlet 
velocity profile encourages the alrstream to separate from the diffuser wall 
or Inner cone. Induces flame stabilization In the separated flow region, and 
causes augmentor components to burn or to Induce combustion Instability. To 
obtain a desirable fuel-air distribution. It Is necessary to understand the 
velocity and temperature profiles In the augmentor. In other words, since the 
velocity and temperature profiles at the exit of the diffuser must be known 
some experimental research was carried out on the annular diffuser of the tur- 
bofan augmentor. 

The test model and the Instrumentation locations are shown In figure 1. 

The test was run at the following conditions: 

(1) Inlet Mach number, 0.08 to 0.60 

(2) Inlet temperature, 300 to 600 K 

(3) Inlet Reynolds number, >5x1 0 4 

(4) Diffuser outlet to Inlet area ratio, 1.54 

(5) Inner to outer annular area ratio, »1.0 

(6) Inner cone, three configurations 

(7) Cone angles, 20°, 30°, and 40° 

The static-pressure recovery coefficients, velocity profiles, and temper- 
ature profiles along the axial distance were measured. The flow separation In 
the diffuser was studied by tuft observation method for low- temperature airflow 
under different Inlet outer-inner velocity ratios and different Inner-cone con- 
figurations. Typical Inlet velocity and temperature profiles of the diffuser 
are shown In figures 2 and 3. 

We found the following results: 

(1) Within the range of the test Mach number (M^ or M e < 0.6), the 
static-pressure recovery coefficients at different axial stations were Indepen- 
dent of Mach number. As shown In figure 4, the Inlet Mach number had little 
Influence on the diffuser Internal flow. Thus, the effect of Inlet Mach number 


c-> 
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can be neglected when testing or analyzing the performance of an annular dif- 
fuser for a turbofan augmentor. 

(2) We also obtained an Important nondlmenslonal flow parameter. 


--(Ve)(vO-(Ve°' 5 )(Vi°- 5 ) 


If we use a diffuser with the same geometry and a similar Inner velocity 
profile and maintain Sq = constant, then the change of static-pressure 
recovery coefficients and nondlmenslonal, combined velocity-temperature combi- 
nation distributions along the diffuser axial distance Is the same, regardless 
of whether the test Is done with low- or high-temperature air, as shown In 
figures 5 and 6. It Is possible to use the low-temperature test to simulate 
real turbofan augmentor diffuser flow; low- temperature airflow separation can 
also be used to understand high- temperature airflow separation. 


RESEARCH ON COMBUSTOR IGNITION PERFORMANCE 

Under the high-altitude relight condition, the combustion chamber Inlet 
pressure can be as low as 0.2 kg/cm2, and the Inlet temperature can be as 
low as -50 °C. Because It Is Important for combustion engineers to understand 
the relation between Ignition performance and combustor Inlet parameters, we 
studied the effect of Inlet parameters on the Ignition performance of combus- 
tors with different fuel-injection systems. 

This test was performed on the high-altitude Ignition test rig shown In 
figure 7. Three combustor configurations were used. The first configuration 
was a gas turbine combustor with a pressure swirl atomizer. (This combustor 
Is similar to a JT-3D engine combustor with one flame tube.) The second con- 
figuration was a one-quarter sector of a full annular combustor with three 
T-shaped vaporizers, as shown In figure 8. The third configuration was similar 
to the second one - the only difference being that the three vaporizers had 
been replaced by three alrblast atomizers. 

A high-energy Ignition plug was used for every configuration. The stored 
energy of the Ignition system was 20 J, 12 J, and 8 J, respectively. The Igni- 
tion procedures were repeated three times under every test condition. If all 
three Ignitions were successful, then each was judged a successful Ignition at 
Its specific condition. 

We drew the following conclusions from this experiment: 

(1) The Inlet pressure has a significant effect on the Ignition perform- 
ance of the combustors. Figure 9 shows the optimum fuel-air ratio as a func 
tlon of Inlet pressure for the three configurations. From figure 9 we can see 
the effect of the Inlet air pressure on Ignition performance of the combustors, 
especially at low pressure. 

(2) The Inlet air temperature has a strong effect on the Ignition perform- 
ance of the combustors. The effects of the Inlet air temperature on maximum 
Ignltable reference velocity and on optimum fuel-air ratio for the three con- 
figurations are shown In figures 10 and 11. As the Inlet air temperature 
decreases, the maximum Ignltable reference velocity decreases. At higher Inlet 
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air temperatures, the effect of the Inlet air temperature on the optimum Ignlt- 
able fuel-air ratio Is very weak; however, at lower Inlet air temperatures the 
effect Is very strong, especially for the vaporizer combustor. 

Figures 12(a), (b), and (c) show the combined effects of Inlet air pres- 
sure and temperature on Ignition performance for the three configurations. 


EXPERIMENTAL STUDY ON AIRFLOW DISTRIBUTION OF THE COMBUSTOR 

The airflow distribution characteristics of a combustor directly affect 
the performance of the combustor such as combustion efficiency, range of sta- 
ble operation, and temperature profile at the combustor exit. A 90° segment 
of the annular combustor was used for this study. The test assembly Is shown 
In figure 13. 

The goals of the study were as follows: 

(1) To determine the effect of Mach number at the combustor Inlet on the 
airflow distribution at an air temperature of 291 K 

(2) To compare the changes of airflow distributions with and without 
burning In the combustor 

Figure 13 shows the eight Instrumentation sections on the test assembly. 
The flow rate through the holes on the wall of the flame tube can be 

calculated as follows: 


G hi - Vhi V 2 ^ 4P 


where 

C . . airflow discharge coefficient of the holes 
dl 

A h1 area of the holes 
Y specific gravity of air 

* 

AP pressure drop through the holes, = P fl ^ - 
g acceleration of gravity 

The flow rate (In percent) through each row of holes that Is relative to the 
total Is 



The relative area (In percent) of the holes of each row of holes Is 
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The following conclusions were obtained from this experimental study: 

(1) At the range tested, there Is no apparent Influence on the airflow 
distribution by the Mach number at the combustor Inlet, as shown In figures 
14(a) and (b). 

(2) Under the conditions 


T 1n ' 291 K 
* / * 

T n . = 1 ~ 2.44 

ex/ in 

the airflow distributions have no obvious change, regardless of whether burning 
Is present or not, as shown In figures 15(a) and (b). However, we are not 
certain that this would be true under real operating conditions. 


EXPERIMENT AL STUDY ON FUEL INJECTORS 

The fuel Injector Is one of the most Important components of a combustor 
or afterburner. In order to determine the effects of geometric parameters, 
airflow parameters, and fuel rate on the fuel-spray characteristics, we have 
made an Initial study of a vaporizer and an alrblast atomizer. 

The test rig assembly Is shown In figure 16. The airflow rate G a 
through the fuel Injector was measured by a turbine flowmeter. The air veloc- 
ity at the Injector exit was 


V a = G a (R)T 2 /AP 

The fuel rate Gf was also obtained by a turbine flowmeter where the air- 
fuel ratio Is given by 


AFR = G a /Gf 

The Sauter mean diameter (SMD) of the fuel spray drops was measured by a 
Malvern laser particle sizer. 


EXPERIMENTAL STUDY ON VAPORIZER FUEL INJECTOR 

Since we tried to understand the spray characteristics - not the vaporiz- 
ing condition - for different vaporizer fuel injectors, we studied the spray 
characteristics under low Inlet air temperature only and did not heat the 
vaporizer. Some of the test vaporizers are shown In figures 17 and 18. 

The following conclusions were made: 

(T) The air velocity V a at the vaporizer exit Is the most Important 
factor affecting fuel atomization. As shown In figure 19, SMD decreases with 
Increasing V a . 
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(2) The effect of air-fuel ratio on fuel atomization Is related to the 
air velocity at the vaporizer exit. If AFR <4 and the air velocity Is low, 
then, as AFR Increases, SMD will decrease. When the air velocity Is much 
higher, the effect Is not apparent, as shown in figures 20 and 21. 

(3) The method and the pressure of injecting fuel Into the vaporizer has 
no significant Influence on fuel atomization, as shown In figure 22. 


EXPERIMENTAL STUDY OF THE AIR-BLAST ATOMIZER 

The geometric parameters of the alrblast atomizers we have tested are 
listed in the following table: 


Number 

Inner channel 
feature 

Outer channel 
feature 

Total exit 
area, 
mm 

0 

Straight flow 

Tangent flow 
holes 

295.4 

11 

With swirl vane 

Tangent flow 
holes 

295.4 

21 

Straight flow 

Rectangular 
flow holes 

359.5 

7 

Straight flow 

Rectangular 
flow holes 

295.4 

17 

Straight flow 

Rectangular 
flow holes 

235.1 


With swirl vane 

Rectangular 
flow holes 

295.4 


lhe test atomizers are shown in figures 23 and 24. 

lhe conclusions were as follows: 

(1) Under the test conditions, SMD decreased with Increasing alrspray 
velocity, as shown In figure 25. 

(2) When AFR <4, SMD decreased with increasing AFR. When AFR >4, there 
was no significant Influence on SMD, as shown in figure 26. 

(3) For the same conditions, SMD decreased as the exit area of the atom- 
izer increased (fig. 27). The reason for such a result has not been determined 
yet. 


There are many interesting subjects related to combustion such as combus- 
tor exit temperature profile, the procedure to control exhaust emission, com- 
bustor diagnostics, swirl combustion, and cooling techniques. We are prepared 
to do further experimental study. 
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Figure 1. - Test section and measurement locations. 
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Figure 2. - Typical, nondlmenslonal diffuser Inlet Figure 3. - Typical, nondlmenslonal 

velocity profile. diffuser Inlet velocity profile. 
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Figure 6. - Nondlmenslonal velocity-temperature comblna 
tlon distribution under constant Sg value for low- 
temperature and high temperature tests. 
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Figure 10. - Effect of inlet temperature 
reference velocity. 
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(b) Experimental combustor with vaporizer. 


Figure 12. - Correlation of maximum ignltable 
reference velocity with inlet parameters. 
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(a) Outer annulus. 
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(b) Inner annulus. 

Figure 15. - Effect of combustion on air-flow 
distribution of annulus. 
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Figure 16. - Test rig assembly. 
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THERMODYNAMICS AND COMBUSTION MODELING 


Frank J. Zeleznik 
Lewis Research Center 
Cleveland, Ohio 


Modeling fluid phase phenomena blends the conservation equations of con- 
tinuum mechanics with the property equations of thermodynamics. The thermo- 
dynamic contribution becomes especially Important when the phenomena Involve 
chemical reactions as they do In combustion systems. The successful study of 
combustion processes requires (1) the availability of accurate thermodynamic 
properties for both the reactants and the products of reaction and (2) the 
computational capabilities to use the properties. A discussion Is given of 
some aspects of the problem of estimating accurate thermodynamic properties 
both for reactants and products of reaction. Also, some examples of the use 
of thermodynamic properties for modeling chemically reacting systems are pre- 
sented. These examples Include one-dimensional flow systems and the Internal 
combustion engine. 


INTRODUCTION 


Modeling fluid phase phenomena always requires a combination of the con- 
servation equations of continuum mechanics with the property functions of 
thermodynamics. For some applications It Is possible to Idealize the model so 
drastically that the continuum equations virtually disappear and the compu- 
tation becomes wholly thermodynamic. In other situations the thermodynamic 
properties are so highly idealized that the computation deals solely with the 
continuum equations. Consequently, thermodynamics and continuum mechanics are 
often viewed as Independent disciplines. However, the complementary nature of 
the two disciplines Is manifested In attempts to generate realistic models of 
physical systems In which combustion occurs. 

The roles of thermodynamics and the conservation equations are easily dis- 
cerned by examining the relevant equations. The conservation equations 
(refs. 1 and 2) for mass, species, energy, and momentum are shown In 
equation (1), where the summation convention Is used for an Index repeated as 
a subscript and a superscript: 

am _ , k. 6m „ k 

it * V mv > ' 0 ‘ IT - m V 


* V m V - v - \ 


a(mu) _ , k. _ # k kj_ 

it * v muu > ■ - V" * ■* v * t Vk 


9Q k. 6U 

m it - v f k ■ m it 


(i) 


v* piv k 1 n _ r k kj 

m tt = r = - mg v.n + t + v.t 


6v_ r k 
6t 


. - pg 1 ^ ♦ a T « 




sk njxsn mr fsi^o 
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In these equations, m Is the mass density, u the Internal energy per unit 

k k 

mass, n the moles per unit mass of species x, v the velocity, q the 

* \ k 

heat flux, y the chemical potential per mole of species x, the diffu- 
sive flux of species x, and R x the volumetric rate of production of species X. 
In the momentum equation, gkl is the metric tensor for three-dimensional 
space, rkj is the stress tensor, fi Is the potential energy, and fk represents 
all volumetric forces over and above those produced by the stress and the 
potential energy. The stress tensor Is composed of a contribution from the 
pressure p and the dissipative stresses. The symbol t represents time, 
v k Is the covariant derivative with respect to the spatial coordinates, and 
6/6t = 3/at +■ vkVfc Is known as the absolute, substantial, or convective 
derivative. The thermodynamic property functions for pressure, Internal 
energy, and the chemical potentials are shown In equation (2). These prop- 
erties may be regarded as functions of mass density, temperature T, and moles 
per unit mass of the species: 


p - p(m,T,n x ) 

u = u(m,T,n x ) (2) 

yX - yX(m,l',n T ) 

Substitution of the thermodynamic Internal energy function Into the con- 
servation equation for Internal energy converts It Into an equation for the 
temperature, and the resulting equation now contains explicit contributions 
from changes In composition and mass density. In reacting fluids the temper- 
ature largely reflects the energy released or absorbed by chemical reactions. 
Thus, we see that the thermodynamic property of Internal energy determines the 
temperature and, consequently, the amount of energy released or absorbed. Of 
course, the temperature, In turn, strongly affects reaction rates and other 
physical properties of the reacting fluid. The thermodynamic pressure equa- 
tion (an equation of state) appears In the equation for the conservation of 
momentum and thus helps to define the local flow field. Clearly, the success- 
ful modeling of combustion processes requires two Ingredients: (1) accurate 

thermodynamic properties for both the reactants and the products of reaction 
and (2) computational capabilities to use the properties In modeling combus- 
tion systems. 

I shall briefly describe some aspects of our work on the two topics of 
property estimation and property utilization. First, I shall discuss the 
estimation of free radical properties In the Ideal gas state and the properties 
of nonideal solutions. Following the discussion of properties, I shall Illus- 
trate some of our computational capabilities to use these properties to model 
the physical systems In which reactions take place. 


ESTIMATION Of FREE RADICAL PROPERTIES 

Free radicals are highly reactive species which are usually present at low 
concentrations In chemically reacting systems. These species play an Important 
role In the chemical reaction mechanisms of combustion. They are Involved In 
both the endothermic (heat absorbing) Initiation steps of the mechanism and the 
exothermic (heat releasing) recombination steps. Recently Bauer and Zhang 
(ref. 3) and others have conjectured that free radicals are also Important 
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species In the mechanism of soot formation In hydrocarbon flames. Hence, the 
thermodynamic properties of free radicals are of considerable Interest. 
Unfortunately, relatively little Is known about the thermodynamic properties 
of radicals over the extended temperature range encountered In combustion 
systems. Hence, we have found It necessary to estimate the thermodynamic 
properties of the phenyl, phenoxy, and biphenyl free radicals which are likely 
to be Important species In the combustion of hydrocarbons and In soot forma- 
tion. Because of their low concentrations and because most combustion processes 
take place at relatively low pressures and high temperatures, It Is usually 
adequate to estimate free radical properties In the Ideal gas state. 

Thermodynamic properties for any species In an Ideal gas state can be 
estimated from the canonical partition function using the formulas of statis- 
tical thermodynamics. This requires knowledge of the species' molecular 
structure, Its vibrational and rotational energy levels, any barriers to 
Internal rotation, and the standard state heat of formation. The difficulty 
Is that this Information Is seldom available for free radicals. As a result 
It becomes necessary to estimate these data from similar, stable species whose 
molecular properties have been measured. These stable, related species are 
called the parent species of the radicals. The species benzene, phenol and 
biphenyl are suitable parent species for the radicals phenyl, phenoxy and 
o-blphenyl, respectively. The relation between these parent species and the 
radicals Is shown In figure 1. This estimation procedure Is based on the work 
of Forgeteg and Berces (ref. 4). In this method It Is usual to assume that the 
bond lengths and bond angles of the radical are Identical to those of the par- 
ent. Also the radical's vibrational frequencies are taken to be the same as 
those of the parent molecule except that three frequencies are eliminated for 
each atom which appears In the parent molecule but not In the radical. In some 
cases, where additional Information about the radical Is available, the assumed 
structure and vibrational frequencies can be adjusted to take Into account this 
extra Information. 

One of the difficulties with this method of estimation Is the selection of 
the appropriate frequencies for elimination. The reason for this Is that a 
molecule's frequencies correspond to vibrations of the molecule as a whole 
rather than vibrations of a particular bond. Naturally this introduces some 
uncertainty Into the resulting properties. To ascertain the degree of uncer- 
tainty one might try eliminating different frequencies or, perhaps, perform the 
estimation with a different parent species. 

Comparisons of the results of our estimation of phenyl and phenoxy prop- 
erties (ref. 5) with estimates made by Benson, et al. (refs. 6 and 7) are shown 
In tables I and II for entropy and heat capacity. Our estimation of the phenyl 
properties from benzene and phenol agree surprisingly well with each other and 
with the room temperature values estimated by Benson. Our estimates of the 
phenoxy properties estimated from phenol also agree very well with Benson's 
estimates which, for this species, were available over an extended temperature 
range. 


NONIOEAL SOLUTION MODELS 

The properties of fuels and oxidizers must be accurately known If com- 
bustion processes are to be modeled properly. Since reactants, especially 
fuels, are often multicomponent, nonideal, liquid solutions, It is desirable 
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to have accurate mathematical representations of the composition dependence of 
their thermodynamic properties. Such representations (solution models) can be 
used to analyze, Interpolate and extrapolate thermodynamic measurements of 
nonideal solutions. Solution models can also be used for the calculation of 
equilibria. The utility of a solution model Is directly proportional to Its 
ability to handle both a wide range of nonideal behaviors and a great diversity 
of multicomponent systems. 

Many solution models have been proposed, but most of them have a limited 
capacity to treat nonideality. Some of them either are only able to deal with 
relatively small deviations from Ideality or cannot handle the complete com- 
position range. Others are limited to binary or ternary solutions. Still 
others exhibit undesirable mathematical behavior. Many contain nonlinear 
parameters which greatly complicate the analysis of experimental measurements. 

I have proposed a solution model (ref. 8) which overcomes most of these dif- 
ficulties. The model Is applicable to solutions containing an arbitrary number 
of constituents. It can treat highly nonideal solutions, electrolyte and non- 
electrolyte, over the complete composition range. It has a relatively simple 
mathematical form, and all parameters appear linearly In the model. Finally, 

It contains both Ideal and regular solutions as special cases. The model's 
expression for the Gibbs free energy Is given In equation (31 where G Is the 
Gibbs free energy, x^ Is the mole fraction of species 1, <p' 1 ) represents 

the elementary symmetric functions, and and are functions of tem- 
perature and pressure. J J 
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I shall Illustrate the flexibility and the capabilities of the model by 
demonstrating how It handles three very different, highly nonideal solutions 
(ref. 9). 


The first solution Is the ethanol-heptane binary system, for which there 
are three sets of excess enthalpy measurements at 30 °C. The solution Is 
endothermic, and measurements extend Into the very dilute solution range. The 
representation of the measurements by the solution model Is shown In figure 2. 
Clearly the model reproduces the data very well over the complete composition 
range. Including the steep gradient that occurs In the dilute ethanol region. 
The second solution Is the chloroform-ethanol binary solution at 50 °C, whose 
excess enthalpy exhibits exothermic behavior over a part of the composition 
range and endothermic behavior over the rest of the range. The representation 

and again the data are well represented by 


A r\ £ A mnr A *) 
III I lyui C J , 


the model. These two solutions display very different behavior yet both can 
be accommodated easily by the same solution model. 


The first two examples were nonelectrolyte solutions. The last example Is 
an electrolyte solution of sodium chloride In water. Measurements for this 
system extend from the very dilute solutions to saturated solutions and cover 
a temperature range of almost 200 °C. A comparison between the measurements 
and the model are shown In figures 4 to 7, where the composition Is given as 
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molality (g mol of solute/kg solvent). Figures 4 to 6 show the Integral heats 
of solution while the excess heat capacity Is shown In figure 7. Once again 
the model gives more than an adequate representation of the experimental 
measurements . 


CHEMICAL REACTION EQUILIBRIA 

Sometimes the behavior of chemically reacting systems can be predicted 
quite well on the basis of chemical equilibrium. Over the years we have 
developed a capability to perform computations reliably for several applica- 
tions Involving chemical reactions (refs. 10 to 16). The calculations are 
based on free energy minimization and assume that the gaseous phase behaves 
Ideally, that the condensed phases are pure liquids and solids, and that the 
volume of the condensed phases Is negligible. The computations can evaluate 
composition and, also, thermodynamic and transport properties for complex 
chemical systems Involving up to 20 chemical elements and 600 reaction 
products . 

Composition and properties can be calculated for a thermodynamic state 
defined by pressure and one of the three variables: temperature, enthalpy, or 

entropy. Alternatively, the thermodynamic state can be specified by assigning 
the volume and one of the three variables temperature. Internal energy, or 
entropy. In addition to the calculation of compositions and properties for an 
assigned thermodynamic state, we can also compute the results for three dif- 
ferent one-dimensional flow processes neglecting heat transfer effects: 

(1) rocket nozzle expansions, (2) Incident and reflected shocks, and 
(3) Chapman-Jouguet detonations. The rocket nozzle and shock calculations can 
be performed for either equilibrium or frozen compositions. Table III Is an 
example of the calculation of the composition and of the thermodynamic and 
transport properties for a simple system composed of hydrogen and oxygen at 
assigned pressures and temperatures. Table IV Illustrates the results for the 
calculation of a one-dimensional shock In a hydrogen-oxygen mixture diluted 
with argon. 


MODELING THE SPARK- IGNITED INTERNAL COMBUSTION ENGINE 

The spark- Ignited Internal combustion engine Is an example of a chemically 
reacting system for which purely thermodynamic predictions are Inadequate. A 
realistic model must take Into account not only thermodynamic properties but 
also heat transfer rates, finite rate chemistry (chemical kinetics), and 
Intake and exhaust flow rates. These effects lead to a complex model and dif- 
ficult calculations. 

A schematic representation of the operation of the spark-ignited Internal 
combustion engine Is shown In figure 8. Each cycle of the operation can be 
described qualitatively by dividing the cycle Into four parts, each part 
spanning ir radians (180°) of crankangle. Each part corresponds approximately 
to one of four processes taking place In the engine. The parts are given names 
which are reasonably descriptive of the four processes: (1) the Intake stroke, 

(2) the compression stroke, (3) the power stroke and (4) the exhaust stroke. 

A fuel-air mixture enters the cylinder during the Intake stroke and Is com- 
pressed by the piston during the compression stroke. The compressed working 
fluid Is Ignited and expanded during the power stroke and Is expelled from the 
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cylinder during the exhaust stroke. Because matter enters and leaves the 
cylinder, the Internal combustion engine Is an open system. Furthermore, the 
burned gases are not completely expelled from the cylinder, and thus the 
residual gases provide a "memory" of each cycle for the subsequent cycle. Each 
cycle occurs rapidly with a cycle repetition rate which ranges from about 25 
to 250 msec. This cycle repetition rate Is comparable to the time scale for 
the complex chemistry which takes place during combustion, and thus the rates 
of chemical reactions become Important. The engine has a movable boundary 
(the piston) and a complex geometry. Finally, the operation of the engine Is 
nonrepeatable, and the engine experiences "cycle-to-cycle" variations. These 
features of the Internal combustion engine make Its modeling difficult Indeed. 
Additionally, often there Is Insufficient Information to describe completely 
the details of the engine operation even If It were mathematically and numer- 
ically possible. 

I have chosen to develop (ref. 17) a hierarchy of models of differing com- 
plexity In order to accommodate the varying amounts of Information available 
In particular cases. The models range from a simple model requiring only 
thermodynamic properties (modeling level 1) to a complex model demanding full 
combustion kinetics, transport properties and poppet valve flow characteristics 
(modeling level 5). The members of the hierarchy are classified according to 
their treatment of several Important features of the Internal combustion engine 
as shown In table V. The description of the operation of the Internal com- 
bustion engine as being divided Into four parts Is adequate for a qualitative 
description. It Is, however. Inadequate to describe the operation of a real 
engine corresponding to modeling at levels 4 and 5. The timing for the Initi- 
ation and cessation of Important events In a more realistic cycle Is sketched 
In figure 9. This figure, often referred to as an Indicator diagram, Is a plot 
of working fluid pressure as a function of cylinder volume for one complete 
cycle spanning 4ir radians (720°) of crankangle. It should be noted that 
between stations 7 and 2 the Intake and the exhaust valves are simultaneously 
open; this corresponds to what Is called the valve-overlap portion of the 
cycle. 


I shall use the fuel-rich combustion of gaseous propane with air as an 
example of modeling the Internal combustion engine at level 5. The air Is 
humidified (75 percent relative humidity), and recirculated exhaust gas Is 
10 percent of the charge. The calculation was a multicycle calculation, but 
only cycle 11 Is shown. The modeling parameters and a summary of results for 
cycle 11 are given In table VI. The combustion of propane was given a full 
kinetic treatment using a mechanism of 121 reactions. Figures 10 to 15 are 
plots of some of the calculated quantities. Straight line segments are appar- 
ent on some of these plots. These are not Indicative of the accuracy of the 
calculations but rather reflect both the number of points saved for plotting 
and the scale of the plots. Figure 10 shows the Indicator diagram for 
cycle 11. Figure 11 shows the fraction of working fluid which has been 
converted to burned, but still reacting, gas during the combustion phase of 


the cycle. The burned and unburned gas 
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(both temperatures were calculated during the combustion phase of the cycle). 
The quenching effect of the endothermic Initiation reactions Is clearly 
visible. The carbon monoxide concentration In mole percent Is shown In figure 
13. Here we can easily see the freezing of carbon monoxide early In the power 
stroke. A similar plot for nitric oxide, as parts per million In mole 
fraction, also shows freezing In figure 14. Finally, In figure 15 Is a 
history of the mass content of the cylinder during cycle 11. 
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TABLE I. - COMPARISON OF PROPERTIES FOR 
PHENYL RADICAL 3 


Property 

Benson 
(ref. 6) 

Estimated 
from benzene 

Estimated 
from phenol 

cO 

s 300 

69.4 

69.03 

69.30 

pO 

C p 300 

18.8 

18.96 

19.51 

r 0 

c p 500 

— 

31.03 

31.70 

r 0 

C p 1000 

— 

46.59 

46.93 

pO 

c p 3000 

— 

59.04 

59.09 


CD-85-17235 

a AII values as eu (cal/g-mol K). 


TABLE II. - COMPARISON OF PROPERTIES 
FOR PHENOXY RADICAL 3 


Property 

Benson 
(ref. 6) 

From phenol 

cO 

5 300 

73.7 

73.70 

rO 

C p 300 

ZZ 5 

22.64 

pO 

C p 400 

29.8 

29.82 

pO 

c p 500 

35.8 

35.80 

pO 

c p 600 

40.6 

40.57 

pO 

c p 800 

47.5 

47.48 

po 

c p 1000 

52.3 

52.17 


a AII values as eu (cal/g-mol K). 
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TABLE III. 

THERMODYNAMIC EQUILIBRIUM PROPERTIES AT ASSIGNED 
TEMPERATURE AND PRESSURE 

CASE NO. 6668 

CHEMICAL FORMULA 
FUEL H 2.00000 
OXIDANT O 2.00000 

0/F= 7.9370 PERCENT FUEL = 11.1894 EQUIVALENCE RATIO= 1 . 0000 PHI= 1.0000 


WT FRACTION ENERGY STATE TEMP 

(SEE NOTE) KJ/KG-MOL DEG K 

1.000000 -9012.332 L 20.27 

1.000000 -12978.762 L 90.18 


THERMODYNAMIC PROPERTIES 


P , MPA 

0 .10132 

0 .10132 

0 .10132 

0.10132 

0 .10132 

0 .10132 

0.10132 

0 .10132 

0.10132 

0.10132 

0.10132 

T, DEG K 

5500 . 0 

5000 . 0 

4500.0 

4000 . 0 

3500.0 

3000.0 

2500 . 0 

2000 . 0 

1500.0 

1000 . 0 

500.0 

RHO, KG/CU M 

1 . 3416-2 

1 .4967-2 

1.7468-2 

2.3376-2 

3.8936-2 

6.2488-2 

8.5005-2 

1.0938-1 

1.4635-1 

2.1955-1 

4 .3909-1 

H, KJ/KG 

55513.9 

52669.0 

47251 . 0 

33833.1 

12360.5 

-1457.79 

-6953 . 36 

-9278.31 

-10745.6 

-11979.6 

-13038.5 

U / KJ/KG 

47961 . 3 

45899.1 

41450.5 

29498.5 

9758.19 

-3079.30 

-8145.34 

-10204.6 

-11438.0 

-12441 . 1 

-13269.2 

G, KJ/KG 

-123677.1- 

107511.2 

-91736.2 

-76960 .8 

-64513.3 

-54702.2 

-46385.9 

-38774.8 

-31609.1 

-24892 . 6 

— 1 8 7 6 7 .8 

S, KJ/CKGHK) 

32.5802 

32.0360 

30.8861 

27.6985 

21.9640 

17.7481 

15.7730 

14.7483 

13.9090 

12.9130 

1 1 . 4586 

M, MOL WT 

6 .055 

6 . 141 

6.450 

7.672 

11.182 

15.383 

17.438 

17 . 951 

18.013 

18.015 

18.015 

( DLV/DLP ) T 

-1 .0081 1 

-1.02147 

-1 . 06320 

-1.15917 

-1 .16898 

-1 . 06245 

-1.01196 

-1 .00121 

-1.00004 

- 1.00000 

- 1.00000 

( DL V/DLT )P 

1 . 0831 

1.2419 

1.7935 

3.2778 

3.8474 

2.2710 

1.2997 

1.0382 

1.0015 

1 . 0000 

1.0000 

CP, KJ/(KG)(K) 

4.6981 

7.2046 

16.3371 

38.8059 

39.0979 

17.2068 

6.6083 

3.3997 

2.6424 

2.2907 

1 . 9548 

GAMMA (S) 

1 .5032 

1 . 3668 

1.2355 

1 . 1639 

1.1268 

1 .1106 

1.1235 

1 . 1705 

1 .2124 

1 .2523 

1.3091 

SON VEL, M/SEC 

3369.4 

3041 . 9 

2677.0 

2246.2 

1712.4 

1341 . 9 

1157.2 

1041.3 

916.2 

760 . 2 

549 . 6 

TRANSPORT PROPERTIES (GASES ONLY) 










CONDUCTIVITY 

IN UNITS OF MI L L IWATTS/ ( CM) ( K ) 








VISC,MILLIPOISE 

1.3674 

1.2723 

1 . 1761 

1.0869 

1 .0207 

0 . 94586 

0.83841 

0.70639 

0.55364 

0 . 37747 

0 . 17324 

WITH EQUILIBRIUM REACTIONS 










CP, KJ/CKGHK) 

4.6844 

7.2044 

16 . 3359 

38.8069 

39. 1072 

17.2071 

6 .6077 

3 . 3935 

2.6423 

2.2907 

1 . 9548 

CONDUCTIVITY 

12.6667 

16.4726 

31.4208 

70 . 9281 

79.3619 

34 .160 1 

9.6648 

3 . 2468 

1 .7354 

0.9740 

0.3603 

PRANDTL NUMBER 

0.5057 

0 . 5564 

0.6115 

0.5947 

0 .5030 

0.4764 

0.5732 

0.7383 

0.8430 

0 .8878 

0.9398 

WITH FROZEN REACTIONS 











CP, KJ/(KG)(K) 

3 . 5289 

3.5100 

3 .4743 

3.3857 

3.2373 

3.1124 

2.9951 

2.8393 

2.6140 

2.2907 

1.9548 

CONDUCTIVITY 

10.3699 

9.6043 

8.7109 

7.4819 

5.8067 

4.3630 

3.3449 

2.4942 

1 .7 049 

0.9740 

0.3603 

PRANDTL NUMBER 

0.4653 

0 .4650 

0.4691 

0.4918 

0.5691 

0.6747 

0.7507 

0 .8041 

0 . 8488 

0 .8878 

0.9398 

MOLE FRACTIONS 












H 

0.66077 

0.65095 

0.61751 

0.50212 

O .25268 

0.05763 

0.00517 

0 . 00012 

0.00000 

O .00000 

0 . 0 0000 

H02 

0.00000 

0 . 00 000 

0 .00000 

0.0 0003 

0.00008 

0.00005 

0 .00 001 

0.000 0 0 

0 .00000 

0.0000 0 

0 . 0 0 0 0 0 

H2 

0 .0 0387 

0.01037 

0 .03230 

0.1 0 047 

0 . 18534 

0 .13479 

0.04283 

0 .00581 

0 . 00 020 

0.0 000 0 

0 . 00000 

H20 

0.00002 

0 . 000 18 

0 . 0 021 3 

0.02961 

0.23353 

0.64341 

0.91089 

0 . 98950 

0 . 99968 

1 . 0 0 0 0 0 

1.00000 

o 

0 . 33098 

0 . 3266 3 

0 . 31048 

0.24985 

0 .11877 

0.02421 

0.00182 

0 .00003 

0.00000 

0.00000 

0 . 00 000 

OH 

0.00363 

0 .00 967 

0 . 02972 

0.08948 

0 .15096 

0.09367 

0.02331 

0 .00214 

0 .000 04 

0.0000 0 

0.0 0 0 0 0 

02 

0 .00073 

0.0 0219 

0 .00786 

0.02845 

0 .05865 

0 . 04624 

0.01596 

0 .00238 

0 .00 009 

0.00000 

0 . 00000 


ADDITIONAL PRODUCTS WHICH WERE CONSIDERED BUT WHOSE MOLE FRACTIONS WERE LESS THAN 0.5Q0Q0E-05 FOR ALL ASSIGNED CONDITIONS 


H202 


03 


H20CS) H2Q(L) 


CD-85-17226 
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TABLE IV. 

SHOCK WAVE PARAMETERS ASSUMING 
EQUILIBRIUM COMPOSITION FOR INCIDENT SHOCKED CONDITIONS 


CASE NO. 

1207 

CHEMICAL I 

FUEL H 

2.0000 0 

FUEL O 

2.00000 

FUEL AR 

1 . 00000 


MOLES 

0.050000 

0.050000 

0.900000 


ENERGY STATE 
CAL/MOL 

13.324 G 
12.952 G 
9.191 G 


TEMP 
DEG K 
300.00 
300.00 
300.00 


0/F = 0.0000 PERCENT FUEL= 100.0000 EQUIVALENCE RATI0= 0.5000 PHI= 0,0000 


INITIAL GAS (1) 
MACH NO. 

3.3528 

3.5052 

3.6576 

3.8100 

3.9624 

4.1148 

4.2672 

4.4196 

4.5720 

4.7244 

Ul, M/SEC 

1100.00 

1150.00 

1200.00 

1250.00 

1300.00 

1350.00 

1400.00 

1450.00 

1500.00 

1550.00 

P , ATM 

0.01316 

0 . 02632 

0.02632 

0.02632 

0 . 02632 

0.02632 

0 .02632 

0.02632 

0.02632 

0 . 02632 

T, DEG K 

300.0 

300.0 

300.0 

300 . 0 

300.0 

300 . 0 

300.0 

300.0 

300.0 

300.0 

RHO, G/CC 

2.0126-5 

4.0252-5 

4.0252-5 

4.0252-5 

4.0252-5 

4.0252-5 

4.0252-5 

4.0252-5 

4.0252-5 

4.0252-5 

H, CAl/G 

0.25457 

0.25457 

0.25457 

0.25457 

0.25457 

0.25457 

0.25457 

0.25457 

0.25457 

0.25457 

U, CAL/G 

-15.578 

-15.578 

-15.578 

-15.578 

-15.578 

-15.578 

-15.578 

-15.578 

-15.578 

-15.578 

G, CAL/G 

-371.94 

-360.97 

-360 . 97 

-360.97 

-360 . 97 

-360.97 

-360.97 

-360.97 

-360.97 

-360.97 

S, CAL/CG) (K) 

1.2407 

1.2041 

1.2041 

1.2041 

1.2041 

1.2041 

1.2041 

1.2041 

1.2041 

1.2041 

M, MOL WT 

37.654 

37.654 

37.654 

37.654 

37.654 

37.654 

37.654 

37.654 

37.654 

37.654 

CP, CAL/CGMK) 

0 . 1372 

0.1372 

0.1372 

0.1372 

0.1372 

0.1372 

0.1372 

0.1372 

0.1372 

0.1372 

GAMMA (S) 

1.6249 

1.6249 

1.6249 

1.6249 

1.6249 

1.6249 

1.6249 

1.6249 

1.6249 

1.6249 

SON VEL, M/SEC 

328.1 

328.1 

328.1 

328.1 

328.1 

328.1 

328.1 

328.1 

328.1 

328.1 

SHOCKED GAS (2) 
U2, M/SEC 

— INCIDENT — EQUILIBRIUM 

666.94 602.16 576.32 

560.44 

549.22 

540.35 

532.52 

524.96 

517.29 

509.40 

P, ATM 

0.10778 

0.27659 

0.32363 

0.36873 

0.41405 

0.46053 

0.50878 

0.55917 

0.61191 

0.66707 

T, DEG K 

1528.1 

1693.3 

1817.1 

1932.4 

2044.1 

2153.0 

2258.8 

2360.3 

2456.7 

2547.3 

RHO, G/CC 

3.3195-5 

7.6873-5 

8.3812-5 

8.9778-5 

9.5277-5 

1.0057-4 

1.0582-4 

1.1118-4 

1.1672-4 

1.2248-4 

H, CAL/G 

91.697 

114.97 

132.65 

149.45 

166.17 

183.16 

200.60 

218.59 

237.18 

256.35 

U, CAL/G 

13.067 

27.830 

39.138 

49 . 982 

60.929 

72.259 

84.165 

96.788 

110.22 

124.45 

G, CAL/G 

-1940.61 

-2079.42 

-2225.77 

-2363.02 

-2496.51 

-2627.06 

-2753.88 

-2875.57 

-2990.89 

-3098.89 

S, CAL/(G)(K) 

1 . 3300 

1.2959 

1.2979 

1.3002 

1.3026 

1.3052 

1.3080 

1.3109 

1.3140 

1.3172 

M, MOL WT 

38.619 

38.617 

38.614 

38.608 

38.597 

38.579 

38.551 

38.510 

38.452 

38.378 

( DLV/DLP ) T 

-1 . 0 00 0 1 

-1.00002 

-1.00005 

-1.00011 

-1.00021 

-1.00041 

-1.00073 

-1.00123 

-1.00193 

-1 .00286 

(DLV/DLT )P 

1 . 0003 

1.0008 

1 .0019 

1.0038 

1.0070 

1.0125 

1.0209 

1.0330 

1.0492 

1.0695 

CP, CAL/CG) (K) 

0 . 1403 

0.1420 

0.1443 

0.1478 

0.1532 

0.1613 

0.1728 

0.1883 

0.2078 

0.2308 

GAMMA (S) 

1 .5795 

1.5700 

1.5576 

1.5403 

1.5165 

1.4859 

1.4497 

1.4110 

1.3734 

1.3400 

SON VEL, M/SEC 

720.9 

756.6 

780.6 

800.6 

817.2 

830.3 

840.4 

848.0 

854.1 

859.9 

TRANSPORT PROPERTIES (GASES ONLY) 

CONDUCTIVITY IN UNITS OF MI L L I CALORI ES/ ( CM) ( K ) ( SEC ) 

VISC,MILLIP0ISE 0.72845 0.77927 0.81611 0.84959 0.88129 

0.91160 

0 . 94048 

0.96775 

0.99326 

1.0169 


WITH EQUILIBRIUM REACTIONS 


CP, CAL/CG) (K) 

0 .1396 

0 . 1415 

0.1431 

0.1478 

0 .1522 

0.1613 

0 .1728 

0 . 1882 

0.2078 

0.2308 

CONDUCTIVITY 

0.1536 

0.1678 

0 .1786 

0.1992 

0.2158 

0.2576 

0 . 3087 

0 . 3799 

0.4721 

0.5823 

PRANDTL NUMBER 

0 . 6619 

0.6573 

0.6539 

0.6304 

0.6215 

0.5706 

0.5264 

0.4796 

0.4371 

0.4031 

WITH FROZEN REACTIONS 










CP, CAt/(G)(K) 

0 .1396 

0.1401 

0.1405 

0.1408 

0 . 1411 

0.1413 

0.1416 

0.1418 

0.1420 

0.1422 

CONDUCTIVITY 

0 .1536 

0.1658 

0 .1747 

0.1829 

0.1907 

0.1984 

0.2058 

0.2131 

0.2201 

0.2270 

PRANDTL NUMBER 

0 .6619 

0.6587 

0.6563 

0.6540 

0.6519 

0.6495 

0.6469 

0.6440 

0.6406 

0.6368 

P2/P1 

8.191 

10.511 

12.298 

14.012 

15.734 

17.500 

19.334 

21.248 

23.253 

25.349 

T2/T1 

5.09* 

5.644 

6.057 

6.441 

6.814 

7.177 

7 . 529 

7.868 

8.189 

8.491 

M2/M1 

1.0256 

1.0256 

1.0255 

1.0253 

1.0250 

1.0246 

1.0238 

1.0227 

1.0212 

1.0192 

RH02/RH01 

1.6493 

1 .9098 

2.0822 

2.2304 

2.3670 

2.4984 

2.6290 

2.7621 

2.8997 

3.0428 

V2(U1-U2)M/SEC 

433.06 

547.84 

623.68 

689.56 

75U .78 

809.65 

867.48 

925.04 

982.71 

1040.60 


MOLE FRACTIONS 


AR 

H 

H02 

H2 

H20 

H202 

O 

OH 

02 


9.2306-1 9.2302-1 
1.2073-7 8.7884-7 
4.6101-8 1.5001-7 
2.6458-6 1.1305-5 
5.1241-2 5.1163-2 
6.699-10 2.448 -9 
2.8572-6 1.2578-5 
7.5326-5 2.0713-4 
2.5622-2 2.5587-2 


9.2294-1 9.2279-1 
4.2747-6 1.5554-5 
2.8252-7 4.7303-7 
3.5151-5 8.8613-5 
5.1019-2 5.0770-2 
4.449 -9 7.237 -9 
3.9852-5 1.0208-4 
4.3641-4 8.0004-4 
2.5527-2 2.5430-2 


4.7175-5 1.2409-4 
7.3544-7 1.0786-6 
1.9569-4 3.8922-4 
5*. 0 3 6 1 - 2 4.9725-2 
1.097 -8 1.574 -8 
2.2885-4 4.6242-4 
1.3432-3 2.1036-3 
2.5285-2 2.5082-2 


9.2144-1 9.2045-1 
2.8829-4 5.9812-4 
1.5029-6 1.9983-6 
7.0587-4 1.1754-3 
4.8791-2 4.7502-2 
2.147 -8 2.796 -8 
8.5400-4 1.4542-3 
3.0994-3 4.3177-3 
2.4818-2 2.4497-2 


9.1909-1 9.1732-1 
1.1195-3 1.9103-3 
2.5453-6 3.1190-6 
1.8088-3 2.5906-3 
4.5833-2 4.3804-2 
3.485 -8 4.173 -8 
2.3021-3 3.4157-3 
5.7144-3 7.2211-3 
2.4133-2 2.3738-2 


ADDITIONAL PKUUUGIS WHICH WbKb tONSlUtKtU BUI WHUSb MULE FKACIIUNS WbRE LESS THAN G.5GGGuE-0o FOR ALL 
ASSIGNED CONDITIONS 

03 H20 ( S ) H20 < L ) 
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TABLE V. - MODELING LEVELS 


Feature 


Modeling level 

5 

4 

3 

2 

1 

Flow 

Poppet valve formulas 

Zero-pressure-drop flow; 

flow reversal only at 0 = 0 or 0 = 37r, 

where discontinuities are possible 

Flame 

Kinetic 

Equilibrium 

Flame 

propagation 

Mass burning formulas 

Instantaneous 

Burned -gas 
chemistry 

Finite rate 

Equilibrium 

Heat 

transfer 

Nonzero 

Zero 


CD-85-17Z32 
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TABLE VI. 


INTERNAL COMBUSTION ENGINE MODEL ZMOTTO 

REF= ZELEZNIK, FRANK J.; AND MCBRIDE , BONNIE J.: 


CYCLE II LEVEL 5 CASE NO. 
MODELING THE INTERNAL COMBUSTION ENGINE. NASA RP-1094, 


111 
1 985 . 


RPM = 3500 0 EGR = 0.100 T(EGR) = 900.6 K SPARK ADVANCE = 20.00 DEG 

MANIFOLD PRESSURE = 0.42760 ATM EXHAUST PRESSURE = 1.00000 ATM 


COMPRESSION RATIO = 10.5 

FUEL PRESSURE = 1.00000 ATM 

UORE S = A 10.9!! N CM OPIC STROKE = 10.312 CM ROD = 17.195 CM CHAMBER AREA =129.030 S9 CM 
TOTAL VOLUME = 1067.33 CC DISPLACEMENT VOLUME = 966.13 CC 


WALL TEMP = 360.0 K 


IVOPEN = 699.00 DEG IVSHUT = 265.00 DEG 

HEAT TRANSFER PARAMETERS Cl = 0.000000 C2 = 0.337200E 00 

KINETIC FLAME FINITE BURNING INTERVAL = 38.7 DEG WIEBE 


EVQPEN = 465 . 00 DEG EVSHUT = 45.00 DEG 

C3 = 0.000000 A = 0.4000 B = 0.8000 

COMBUSTION TAU =0.8930E-04 SEC BETA = 3.200 

WT FRACTION ENERGY STATE 


CHEMICAL FORMULA 
FUEL C 3.00000 H 8.00000 

AIR N 1.56168 0 0.41959 

* AIR INCLUDES 

A/F= 12.5823 


MASS PER CYCLE (G) 


AR 0.00936 C 0.00032 
0.02348 MOLE FRACTION WATER 


ENERGY 
CAL/MOL 
1.000000 -24821.770 

*1. 000000 -1386.340 

RELATIVE HUMIDITY= 0.7500 


TEMP 
DEG K 
298.15 
298.15 


PERCENT FUEL= 7.3625 EQUIVALENCE RATIO= 1.2500 
PERFORMANCE PARAMETERS FOR ONE CYLINDER 
MEAN INLET MASS FLOW RATE CG/SEC) 


PHI= 1.2647 


MEAN EXHAUST MASS FLOW RATE (G/SEC) 


TOTAL 

FUEL 

AIR 


0.49685 
0 . 02860 
0.35987 


CHARGE 

FUEL 

AIR 

NET 


12.6027 

0.8351 

10.5073 

12.6376 


EXHAUST 

CO 

NOX 

NET 


13.2527 
0.71753 
0 .0001 3 
13.2526 


ENERGY PER CYCLE (JOULES) 

INDICATED WORK 404.672 
INDICATED PUMP WORK -86.886 
HEAT LOSS 352.020 
CHEM. ENERGY 1305.823 


AVERAGE ENERGY RATE - POWER (KW) 

INDICATED POWER 11.803 

INDICATED PUMP POWER -2.534 

HEAT LOSS RATE 10.267 

EXHAUST POWER 16.593 

MISCELLANEOUS 


CYCLE EFFICIENCIES 


NET WORK 
HEAT LOSS 
EXHAUST 


0.243361 

0.269577 

0.435675 


INDICATED MEAN EFFECTIVE PRESSURE (ATM) 
PUMP MEAN EFFECTIVE PRESSURE (ATM) 

MEAN TORQUE (NEWTON-METERS) 

COMPOSITE EXHAUST GAS MOLE FRACTIONS AT 900.63 K AND I. 

AR 0.007952 CH4 0.000914 CO 

C2H4 0.000131 H2 0.026569 H20 

02 0.000070 

FRESH CHARGE MOLE FRACTIONS AT 292.50 K AND 0.4276 ATM 
AR 0.000861 CH4 0.000099 CO 

C2H4 0.000014 H2 0.002878 H20 

C3H8 0.043867 AIR 0.847803 

NOTE: INLET AND EXHAUST VALUES CALCULATED WHEN THE VALVES CLOSE. 

COMPUTER CYCLE TIME- 48.858 SEC 


4.1338 
-0.8876 
25.2886 

MOLECULAR WEIGHT = 26.947 


0.052088 

C02 

0.077497 

C2H2 

0.002239 

0.169476 

NO 

0.000009 

N2 

0.663041 

MOLECULAR 

WEIGHT = 

29.192 



0.005643 

C02 

0.008395 

C2H2 

0.000243 

0.018359 

N2 

0.071827 

02 

0.000008 
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PHENYL 
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PHENOXY 
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O-BIPHENYL 
H H H 

H HH H 
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H 



H 


BIPHENYL 
H H H H 
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Figure 1. - Radicals and parent species. 



Figure 2. - Excess enthalpy of ethanol (1) - heptane 
(2) binary system at 30 °C. 
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INTEGRAL HEAT OF SOLUTION, AH jnt , cal/g-mol 



Figure 3. - Excess enthalpy cf chloroform (1) - 
ethanol (2) binary system at 50 °C. 



Figure 4. - Integral heat of solution of sodium 
chloride in water. 
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INTEGRAL HEAT OF SOLUTION, AH jnt , cal/g-mol INTEGRAL HEAT OF SOLUTION. AH jnt , cal/g-mol 



Figure 5. - Integral heat of solution of sodium 
chloride in water. 



Figure 6. - Integral heat of solution of sodium 
chloride in water. 
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OTTO CYCLE INDICATOR DIAGRAM 
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EFFtCT OF FLAME-TUBE HEAD STRUCTURE ON COMBUSTION CHAMBER PERFORMANCE 


Gu Minqqi 

Shenyang Aero-Engine Research Institute 
Shenyang, People's Republic of China 


This paper presents the experimental combustion performance of a pre- 
mixed, pilot-type flame tube with various head structures. The test study 
covers an extensive area: efficiency of the combustion chamber, quality of 

the outlet temperature field, limit of the fuel-lean blowout, ignition per- 
formance at ground starting, and carbon deposition. As a result of these 
tests, a nozzle was found which fits the premixed pilot flame tube well. The 
use of this nozzle optimized the performance of the combustion chamber. 

The tested models had premixed pilot chambers with two types of air-film- 
cooling structures, six types of venturi-tube structures* and secondary fuel 
nozzles with two small spray-cone angles. 


INTRODUCTION 

After a period of engine operation, a premixed pilot-type flame tube 
(fig. 1) with a pressurized fuel nozzle may have difficulty in starting. 
Sometimes the pilot chamber may not ignite at all; this seriously affects the 
starting performance of the engine. This problem is generally believed to be 
due to carbon deposition in the premixed pilot chamber (8 g carbon deposited 
in a single premixed pilot chamber after 190 min of running). In the case of 
high-load running, the gas temperature in the center of the premixed pilot 
chamber can be raised as high as 1200 °C; this can make the high-density fuel 
spray crack and deposit on the downstream section of the nozzle. In addition, 
the wall temperature would be high, causing the fuel spray to form coke on the 
wall, which (coupled with fuel deposition) would make the wall rough. Then the 
carbon particles would hit the cores of deposition and would be deposited on 
the wall. Thus, the wall becomes uneven and very rough. When the fuel spray 
again hits the wall and flows on its uneven soft carbon layer, the fuel film 
is also very uneven (unlike the film on a smooth metal wall). As a result, the 
spray quality becomes very poor and the engine fails to start. 

The objective of this experiment was to Improve the ground starting per- 
formance while not damaging the other types of combustion performance. The 
study concentrated on the effect of premixed pilot chamber structures and 
nozzle spray angles on combustion performance. 

This paper presents a study of the effect on combustion of (1) pre- 
mixed pilot chambers with ai r-f 1 Im-cooling and venturi-tube structures and 
(2) small-angle secondary nozzles. It emphasizes their effects on ground 
starting performance and on carbon deposition within the pilot chambers. 


TEST MODELS 

Premixed Pilot Chamber With Air-Film-Cooling Structure 

It is known that the ai r-f ilm-cool ing structure enables free carbon 
particles to be carried away by the film-cooling air flow, thus reducing or 
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avoiding the contact of these particles with the wall. This helps reduce the 
wall temperature. Two models were evaluated (fig. 2). Eighty film-cooling 
holes (4*2) were added to the wall of these two premixed pilot chambers. The 
ratio of the cooling-hole opening area to the total opening area in the flame 
tube wall was 0.022. The air velocity in the outlet of the air-film-cooling 
ring was 15 to 26 m/sec. 


Premixed Pilot Chamber With Venturi- Tube Structure 

A venturi tube was mounted in the rear of the swirler (fig. 3). It made 
the fuel injected from the nozzle hit the venturi-tube wall and then mix with 
the air coming from the inner section of the swirler, thus avoiding contact of 
the fuel with the flame tube wall. In addition, because the venturi tube was 
used, the position of the recirculation zone moved further downstream of the 
nozzle; as a result, the temperature of the carbon-generating zone downstream 
of the fuel nozzle was reduced. This helped reduce the carbon particles. 

Six venturi tubes of different sizes were evaluated. (See table I and 
fig. 4.) The ratio of the outer to the inner area of the swirler was 1.16 
when the venturi tube was mounted. 


Secondary Nozzle With Small Spray Angle 

fuel entered from the centrifugal nozzle in the ignition process in order 
to reduce or eliminate the probability of the sprayed fuel hitting the- wall. 
Two kinds of small spray-angle secondary nozzles were studied. (See table II 
and fig. 5.) The spray angle was made small enough that the spray would not 
hit the inner wall of the premixed pilot chamber. Thus, the fuel spray cone 
was not impaired. The fuel-air ratio was appropriate near the Ignitor, and a 
good starting performance was maintained (even if there was carbon deposition 
in the pi lot chamber) . 


RESULTS AND ANALYSIS 

Amount of Carbon Deposition in Premixed Pilot Chamber 

The experimental results for various models are listed in table III. 
(Results for the original flame tube are included for comparison.) The 
results indicate that the venturi tube is the best structure for eliminating 
carbon deposition. There was almost no carbon deposition in the pilot chamber, 
and there was only a small amount of carbon deposition in the 20- to 30-mm 
position downstream of the pilot chamber exit. The objective of reducing car- 
bon deposition in the pilot chamber was achieved. 

Results from the water analog test indicate that there was no recircula- 
tion zone in the venturi tube. The recirculation zone in the premixed pilot 
chamber was located far downstream of the nozzle (fig. 6). The heat radiation 
from the flame-tube center to the nozzle was blocked outside of the venturi 
tube, thus reducing the temperature downstream of the nozzle. This had a 
determinate effect on the reduction of the carbon particles. It was for 
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this reason that the amount of carbon deposition in the flame tube with a ven- 
turi tube was obviously reduced. 

Adding an air-film-cooling structure in the premixed pilot chamber helped 
reduce carbon deposition. The reduction rate was usually as much as 20 to 30 
percent near the air-film exits. But, at the exits of the air films, the air 
velocity was small, so was the air-flow momentum. Nevertheless, the momentum 
of the fuel injected from the nozzle was rather large, so the film flow could 
not prevent the fuel droplets from hitting the wall. It is only possible for 
the film air to bring away free carbon particles. For this reason, this 
structure is not very effective in reducing carbon deposition. 


Ground Starting Performance 

The results for five models at the ground starting ignition condition are 
shown in figures 7 and 8 where V 2 is the inlet velocity of the combustion 
chamber. The curves of figure 7 indicate that there is an optimal value for 
each curve. This velocity ranges from 15 to 20 m/sec. The ignition perform- 
ance was best for the flame tube with a model I venturi tube. 

Figure 8 indicates that the ground starting Ignition performance Improved 
greatly when model I and II nozzles were used under carbon deposition condi- 
tions in the premixed pilot chamber. At 17.5 m/sec of inlet velocity the 
fuel-lean ignition limit Improved; the limiting air-fuel ratio Increased from 
36.8 to 63.2. It is readily seen from these results that, under the carbon 
deposition condition in the premixed pilot chamber, adoption of a small-angle 
secondary nozzle obviously Improves the ground starting ignition performance. 

In addition, it is known from the curves that when carbon is deposited in 
the premixed pilot chamber the air-fuel ratio of the fuel-lean ignition limit 
decreases from 47 to 36. This confirms that carbon deposition in the premixed 
pilot chamber deteriorates the ground starting ignition performance. 


Outlet Temperature Field of Combustion Chamber 

Table IV shows the test results of various venturi-tube models. It is 
readily seen from the table that the size of the venturi tube mounted in the 
premixed pilot chamber has a great effect on the combustor outlet temperature 
field. The outlet temperature distribution factor (OTDF, table IV) of the 
model V combustion chamber was 0.155; this was much better than that for the 
original model (0.229). The OTDF's of models I and II were similar to that of 
the original model, but those of models IV and VI were inferior to that of the 
original model. 

The convergent angle of the front section and the outlet angle of the 
venturi tube have an obvious effect on the combustor outlet temperature dis- 
tribution. When the convergent angle of the front section of the venturi tube 
increases, the uniformity of the combustor outlet temperature field decreases. 
When the convergent angle was Increased from 27°34' to 67°22', the OTDF 
Increased from 0.209 to 0.301. When the venturi-tube outlet angle became 
rather small, excessive fuel concentrated in the center of the flame tube; 
this caused the temperature to rise in the center region and the outlet 
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temperature distribution to deteriorate. For example, reducing the outlet angle 
of the venturi tube from 88° to 56° Increased the OTDF from 0.155 to 0.301. 

Adoption of a secondary nozzle with a small spray angle also Improved the 
uniformity of the combustor outlet temperature field - the OTDF of model II was 
0.175, the radial temperature distribution factor (RTDF, table IV) was 0.077. 


Stability of Combustion Chamber 

Figure 9 shows experimental results of the fuel-lean blowout limit for 
various models. The shapes of the curve are nearly the same for various 
models. The fuel-lean blowout limit for the model II venturi tube was wider 
because this venturi tube had no throat section; thus, a small recirculation 
zone was formed Inside the venturi tube. The size of the recirculation zone 
was limited to the nozzle diameter range (fig. 10). 

The fuel-lean blowout limits of the model III, IV, and V venturi tubes were 
narrower. It Is known from water analog test results that there Is a high- 
speed rotating flow In the throat section of a venturi tube (fig. 11). This 
unstable flow deteriorates the fuel-lean blowout limit. 

When a secondary nozzle with a small spray angle was adopted, the fuel- 
lean blowout limit Improved In the wide velocity range. 


CONCLUSIONS 

This study led to the following conclusions: 

(1) Adoption of the secondary nozzle with a small spray angle not only 
solves the engine's starting Ignition problem when there Is carbon deposition 
In the premixed pilot chamber but also extends the ground starting Ignition 
range and Improves the uniformity of the combustor outlet temperature field. 

(2) Adoption of an alr-fllm-coollng structure reduces the amount of carbon 
deposition on the wall of the premixed pilot chamber. But since there Is only 
limited cooling air available, It Is not easy to totally eliminate the carbon 
deposition In the premixed pilot chamber. 

(3) Mounting the venturi tube In the rear of the swlrler Is effective In 
eliminating carbon deposition In the premixed pilot chamber. But the geometry 
of the venturi tube affects the uniformity of the combustor outlet temperature 
field and the ground starting Ignition performance; therefore, characteristic 
geometric parameters of the venturi tube have to be elaborately adjusted. 
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TABLE I. - CHARACTERISTIC DIMENSIONS OF VENTURI TUBE 


Model 

Inlet 

Diameter, 

mm 

Length, 

mm 

Throat 
diameter , 
mm 

Outlet 

cone 

angle, 

deg 

Outlet 

diameter, 

mm 

I 


36 

30 

56° 

40 

MM 

■■■ 


25 

42 

36°32 ' 

45.3 


K 


25 

30 

56° 

40 


n 


20 

30 

56° 

35 

wm 



20 

30 

88° 

40 

mam 

1 

II 

15 

15 

0° 

15 


TABLE II. - CHARACTERISTICS OF NOZZLE 


Model 

Spray 

angle, 

a 

sec , 
deg 

Fuel 

flow, 

G f 

sec’ 

liter/hr 

“double, 

deg 

G f 

double, 

liter/hr 

Fuel supply pressure, 

P-j =10 kg/cm 2 

sec 

Fuel supply pressure, 

2 

Py = 30 kg/cm 

double 

Original 

91 

24.6 

103 

602.4 

I 

51 

24.6 

100 

602.4 

II 

62 

23.16 

100 

593.4 
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TABLE III. 


EFFECT OF VARIOUS MODELS 


ON CARBON DEROSI I ION 

[Test conditions at inlet of combustion chamber: 
air pressure, P 2 , 6 kg/cm 2 ; air temperature, 
T 2 , 595 K; air flow, G a , 3.135 kg/sec; air- 
fuel ratio, a/f, 70; test duration, 2.5 hr.] 


Model 

Amount of 
carbon 
deposition, 
9 

Relative 

percentage 

Original flame tube 

1 .3 

100 

Model I venturi tube 

.072 

5.54 

Model VI venturi tube 

0 

0 

Model I air film cooling 

1.076 

82.8 

Model 11 air film cooling 

.932 

71 .7 


TABLE IV. - EFFECT OF VENTURI-TUBE GEOMEIRY ON PERFORMANCE OF COMBUSTION CHAMBER 


[lest conditions at inlet of combustion chamber: air pressure, P 2 , 

6 kg/cm2; air temperature, T-?, 695 K; air flow, G a , 3 kg/sec; 
air fuel ratio, a/f, 52.2.] 


Model 

Original 

I 

11 

III 

IV 

V 

VI 

Inlet convergent angle 



27°34 ' 

.. ... „ _ 

62° 

67°22 1 

67°22 1 

122 ° 

Outlet cone angle 

' 

56° 


56° 

56° 

88 ° 

— 

a 0TDF 

0.229 

0.209 

0.223 

0.236 

0.301 

0.155 

0.406 

aRTDF 

.083 

.143 

.098 

.133 

.138 

.085 

.110 

a n c 

.984 

.984 

.975 

.988 

.997 

.979 

.929 


a 0T0F , (15 - T 2 ) /( 1 3 - T£); RTDF = (1 3 - T$)/(T$ - T^) ; and 

max av max(h) av 

n = efficiency of combustion chamber where T* is the average combustion 
c 3 av 

outlet temperature, T* the maximum combustor outlet temperature, and T* 

3 max 3 max(h) 

the maximum combustor outlet temperature along the blade height. 
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(e) Model V. (f) Model VI 

Figure 4. - Concluded. 



Figure 5. - Secondary nozzle. 




Figure 6. - Flow spectrum In flame- tube head 
with model I venturi tube. 



0 5 10 i r ' 20 25 30 V 2 m/soc 


Figure 7. - Effect of venturi-tube size on 
ground Ignition performance. Test condi- 
tions at combustion chamber Inlet: air 

pressure, p£, 1.04 kg/cm?; air tempera- 
ture, T 2 , 331 K. 
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Figure 8. - Effect of fuel spray angle on 
ground Ignition performance with carbon 
deposition In premixed pilot chamber. 
Test conditions at combustion chamber 
Inlet: air pressure, Pji, 104 kg/cm 2 ; 

air temperature, T|, 311 K. 
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Figure 9. - Fuel-lean limits of various models. 
Test conditions at combustion chamber Inlet: 
air pressure, Pp, 4.75 kg/cm?, air tempera- 
ture, T 2 , 655 K. 
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EXPERIMENTS AND MODELING OF DILUTION JET FLOW FIELDS 

James D. Holdeman 
NASA Lewis Research Center 
Cleveland, Ohio 


This paper presents experimental and analytical results of the mixing of 
single, double, and opposed rows of jets with an Isothermal or variable- 
temperature main stream In a straight duct. This study was performed to Inves- 
tigate flow and geometric variations typical of the complex, three-dimensional 
flow field In the dilution zone of gas-turbine-engine combustion chambers. 

The principal results, shown experimentally and analytically, were the 
following: (1) variations In orifice size and spacing can have a significant 

effect on the temperature profiles; (2) similar distributions can be obtained. 
Independent of orifice diameter, If momentum-flux ratio and orifice spacing are 
coupled; (3) a first-order approximation of the mixing of jets with a variable- 
temperature main stream can be obtained by superimposing the main-stream and 
jets-ln-an-lsothermal -crossflow profiles; (4) the penetration of jets Issuing 
from slanted slots Is similar to that of jets from circular holes, but the 
mixing Is slower and Is asymmetric with respect to the jet centerplanes, which 
shift laterally with Increasing downstream distance, (5) double rows of jets 
give temperature distributions similar to those from a single row of equally 
spaced, equal-area circular holes; (6) for opposed rows of jets, with the ori- 
fice centerlines In line, the optimum ratio of orifice spacing to duct height 
Is one-half the optimum value for single-side Injection at the same momentum- 
flux ratio; and (7) for opposed rows of jets, with the orifice centerlines 
staggered, the optimum ratio of orifice spacing to duct height Is twice the 
optimum value for single-side Injection at the same momentum-flux ratio. 

In Illustrating these results, the mean temperature measurements are com- 
pared with profiles calculated using an empirical model based on assumed verti- 
cal profile similarity and superposition and with distributions calculated 
using a three-dimensional elliptic code that had a standard k-e turbulence 
model. The empirical model predictions are very good within the range of the 
generating experiments, and the numerical model results, although they exhibit 
too little mixing, correctly describe the effects of the principal flow and 
geometric variables. 


A. / A 
j m 


C 

C 


SYMBOLS 

jet to main-stream area ratio, A^/A m = ( -ir/4) /[ ( S/H Q ) (H Q /D) 2 ] 
for one-side Injection and ( ir/2) /[ ( S/H Q ) (H Q /D) 2 ] for two-side 
Injection 

< s/ v yfT, eq. (3) 
orifice discharge coefficient 
orifice diameter 
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DR jet to main-stream density ratio, DR s T /T, 

m j 

H q duct height 

2 

J jet to main-stream momentum-flux ratio, J = (DR)R 

M jet to main-stream mass-flux ratio, M = ( DR ) R 

R jet to main-stream velocity ratio, R = V^/U^ 

S spacing between orifice centers 

S x spacing between orifice rows 

7 temperature 

7j jet exit temperature 

T m main- stream temperature 

U velocity 

U main-stream velocity 

m 

Vj jet velocity 

W l/2 half-widths above (+) or below (-) the centerline (ref. 5) 

Wj/w-j jet to total mass-flow ratio, Wj/wj = ^{QR) J (C d ) 

(Aj/An,)/(l + y/W)i (C d ) (Aj/A m )J 

x downstream coordinate (x = 0 at injection plane) 

y cross-stream (radial) coordinate (y = 0 at wall; y = y c at location 

of minimum temperature in a line x = constant and z = constant) 

z lateral (circumferential) coordinate (z = 0 at centerplane) 

e (T m - T)/(T m - Tj), eq. (1) 

e c temperature difference ratio at y c 

4 - 

e~. minimum temperature difference ratio above (+) or below (-) the 
centerline, (fig. 4) 


INTRODUCTION 

The problem of jets in crossflow has been rather extensively treated in 
the literature, to the point that it can almost be called a classical three- 
dimensional flow problem. Although these studies have all contributed to the 
understanding of the general problem, the information obtained in any given 
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study Is naturally determined by the motivating application and therefore may 
not satisfy the specific needs of diverse applications. 

Consideration of dilution zone mixing In gas-turbine combustion chambers 
has motivated several previous studies of the mixing characteristics of jets 
Injected normally Into a ducted crossflow (refs. 1 to 12). One factor making 
the combustor dilution zone jet-ln-crossf low application unique Is that It Is 
a confined mixing problem - 10 to 50 percent of the total flow enters through 
the dilution jets. The result Is that the equilibrium temperature of the exit- 
ing flow may differ significantly from that of the entering main-stream flow. 

To control or tailor the combustor exit temperature pattern, we must be able 
to characterize the exit distribution In terms of the upstream flow and geo- 
metric variables. This requires that the entire flow field be either known or 
modeled . 

Empirical correlation of experimental data can provide an excellent pre- 
dictive capability within the parameter range of the generating experiments 
(e.g., refs. 4 to 6), but empirical models must be used with caution, or not 
at all, outside this range. Physical modeling, in various levels of sophisti- 
cation and complexity, may be used to obviate this weakness. In this regard, 
several one and two-dimensional. Integral and differential jet-ln-crossf low 
models have been developed and shown to give, for example, trajectory predic- 
tions that are in good agreement with experiments. These models may provide 
Insight into the dominant physical mechanism(s) and may predict some of the 
characteristic parameters well, but they rarely provide sufficient Information 
to quantify the flow field In three coordinate directions. 

Recently, rapid advances have been made in the capability of computational 
fluid dynamics models and In their application to complex flows such as jets 
In crossflow (refs. 13 to 16). These models are, however, still In the devel- 
opment and verification stage. They have been shown to be capable of predict- 
ing trends in complex flows, but their capability to provide accurate, 
quantitive, and grid-independent calculations of these flows has not yet been 
demonstrated (refs. 14 to 16). 


PREVIOUS RESULTS AND THE CURRENT STUDY 

The study in references 17 to 19 was performed to extend the available 
experimental data on, and the empirical correlations of, the thermal mixing of 
multiple jets in crossflows so that they would include geometric and flow 
variations characteristic of gas turbine combustion chambers - namely, variable 
temperature main stream, flow area convergence, noncircular orifices, and 
double and opposed rows of jets. These experiments are a direct extension of 
those in reference 1. The effect on the dimensionless temperature distribu- 
tions of varying the jet to main-stream density ratio, the momentum-flux ratio, 
and the orifice size and spacing are presented in reference 2. 

From the data in reference 1, an empirical model was developed (refs. 4 
and 5) for predicting the temperature field downstream of a row of jets mixing 
with a confined crossflow. The effects of separately varying the Independent 
flow and geometric variables and the relations among these variables which 
optimize the mixing are reviewed In reference 12. This study was conducted 
using an interactive microcomputer program that is based on the empirical model 
of reference 6. 
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The resulls of these Investigations may be summarized as follows: (1) 

mixing Improves with Increasing downstream distance; (2) the momentum-flux 
ratio Is the most significant flow variable; (3) the effect of density ratio 
Is small at constant momentum flux ratio; (4) at any given momentum-flux 
ratio, decreasing orifice spacing while maintaining a constant orifice diameter 
reduces penetration and Increases lateral uniformity; (5) Increasing orifice 
diameter while maintaining a constant spacing- to-dlameter ratio Improves pene- 
tration but Increases lateral nonuniformity; (6) Increasing orifice diameter 
at a constant orifice spacing Increases the magnitude of the temperature dif- 
ference, but jet penetration and profile shape remain similar; (7) profiles for 
conditions with equivalent coupling of orifice spacing and momentum flux ratio 
show similar distributions; and (8) smaller momentum- flux ratios (larger spac- 
ing) require a greater downstream distance for equivalent mixing. 

Results from the experiments In references 17 to 19 that are considered 
In this paper Include the effects of variations In orifice size and spacing, 
coupled spacing and momentum- flux ratio, variable temperature main stream, 
noncircular orifices, and double and opposed rows of Inline and staggered jets. 
Also, temperature field measurements from several experiments are compared with 
distributions calculated using an empirical model based on assumed vertical 
profile similarity and superposition (refs. 17 to 19) and using a three- 
dimensional elliptic code with a standard k-e turbulence model (ref. 14). The 
results show the capability of these models to predict the effects of the 
principal flow and geometric variables. 

A more complete presentation of the experimental results and a discussion 
of the empirical modeling performed In this study are given In references 1, 

4, and 17 to 19. Selected experimental and analytical results from these 
studies and from reference 14 are also given In references 2, 5, and 20 to 23. 


FLOW FIFI.D DESCRIPI ION 


Figure 1 shows a schematic of the dilution jet flow field for jet Injec- 
tion from the top wall. The temperature field results are presented In three- 
dimensional oblique views of the temperature difference ratio 6: 


A sequence of experimental profiles of this parameter at several loca- 
tions downstream of the Injection plane Is shown In figure 2. In the three- 
dimensional plots the temperature distribution is shown In the y,z-planes 
normal to the main flow direction x. The coordinates y and z are, respec- 
tively, normal and parallel to the orifice row. Note that the jet fluid is 
identified by the larger values of the e parameter (l.e., e = 1 If T = lj. 
and e = 0 if 1 = T m ) . The equilibrium e for any configuration is equal 
to the fraction of the total flow entering through the dilution jets Wj/Wf. 

The orifice configurations Investigated are shown In figure 3. The pri- 
mary independent geometric variables for each orifice configuration are the 
spacing between adjacent orifices S, the orifice diameter D (for noncircular 
orifices, this is taken as the diameter of a circle of equal area), and, for 
double rows, the axial spacing between rows Sx- These are expressed In 
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dimensionless form as the ratio of the orifice spacing to duct height S/Hq, 
the ratio of the duct height to orifice diameter Ho/D, and the ratio of the 
axial spacing to the duct height S x /Ho. 


EXPERIMENTAL CONSIDERATIONS 

The dilution-jet mixing characteristics were determined by measuring tem- 
perature and pressure distributions with a vertical-rake probe positioned at 
different axial and lateral stations. This probe had 20 thermocouple elements. 
A 20-element total-pressure rake and a 20-elemerit static-pressure rake were 
located nominally 5 mm (0.05 Hq) on each side of the thermocouple rake. The 
center- to- center spacing between sensors on each rake was 0.05 Ho. 

This probe traversed a matrix of 48 to 64 z.x-plane survey locations. 

The flow field was mapped In the z-dlrection over a distance of 1 or 1.5 times 
the hole spacing S at Intervals of S/10. For most tests, the x,y-plane 
containing the orifice centerline (centerplane) was at the center of the span 
surveyed; that Is, data surveys were from midplane to midplane. 

Measurements In the x-dlrectlon were made at up to five planes with 
0.25 < X/Hq < 2. Because the objective In this application is to Identify 
dl lution- zone configurations that will provide a desired mixing pattern within 
a given combustor length, the downstream stations are defined In Intervals of 
the duct height Hq rather than Intervals of the orifice diameter D. 


FLOW FIELD MODELS 


Empirical 


The empirical model for the temperature field downstream of jets mixing 
with a confined crossflow Is based on the observation that properly nondlmen- 
slonallzed vertical temperature profiles can be expressed In the following 
self-similar form (ref. 5) for any location in the flow field: 


e 
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e m1n 

6 mfn 


exp 
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(y - y c ) 


w 
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where © is the temperature difference ratio at vertical location y, and 

— 

e m1n’ e min’ W l/2’ W l/2’ e c ’ an< ^ y c are sca "* in 9 parameters as shown In 

figure 4. Correlations have been developed for each of these In terms of the 
Independent variables J, S/D, Hq/D, Z/S, and x/Hq. The correlations In 
reference 5 for a single row of jets In a uniform temperature crossflow have 
been extended for predicting the temperature field downstream of single, dou- 
ble, or opposed rows of jets, either Inline or staggered. Injected Into an 
Isothermal or nonl sothermal main stream, with or without flow-area convergence 
( refs . 17 to 19). 
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Numerical 


The numerical code used In this Investigation Is based on the USARTL 
three-dimensional, fully elliptic, turbulent flow model (ref. 24) and uses 
pressure and velocity as the main hydrodynamic variables. This code, or simi- 
lar versions thereof, has been used In previous validation and assessment 
studies (refs. 14 to 16). 

The governing equations are represented by finite-difference approxima- 
tions on a staggered grid system. The differencing technique employed Is 
hybrid for convective terms with central differencing of all other terms. The 
velocity-pressure coupling Is handled by the SIMPLE algorithm of Patankar and 
Spalding (refs. 25 and 26). Uniform velocities and mass flow rates were used 
at all Inflow boundaries. The code contained a conventional k-e turbulence 
model, and standard values of the constants Cp, C-| , and C 2 were used 
(l.e., Cp = 0.09, Ci = 1.44, C 2 = 1.92). The rms turbulence Intensity was 
chosen to be 4.5 percent of the local mean velocity, the Inlet length scale was 
2 percent of the jet diameter and duct height for the jet and main stream, 
respectively, and the turbulent Prandtl number was 0.9 for all calculations. 


RESULTS AND DISCUSSION 

The following paragraphs describe the experimental results and compare 
them with those of the empirical and numerical model calculations, and Illus- 
trate the effects of the primary Independent variables. The flow and geometry 
conditions corresponding to the figures shown are given In table I. 


Orifice Size and Spacing 

At constant orifice area, changes In orifice size and spacing can have a 
significant Influence on the 6 profiles. This Is shown by the experimental 
profiles In figure 5 where jets from closely spaced small orifices underpene- 
trate and remain near the Injection wall (fig. 5(a)), and jets from widely 
spaced larger orifices overpenetrate and Impinge on the opposite wall 
(fig. 5(b)). In this figure, and In several others In which the orifice 
spacing Is different for different parts of the figure, a duct cross-section 
Is Included to Indicate the region for which data are shown. 


The data for these conditions, at x/Hp = 0.5, are compared with calcu- 
lated distributions In figure 6. The empirical model reproduces the data very 
well In the small orifice case because the data are consistent with the major 
assumption In the empirical model, namely that all vertical temperature dis- 
tributions can be reduced to similar Gaussian profiles. The empirical model 
does not do as well In the larger orifice case, however, because the Impinge- 
ment of the jets on the opposite wall results in vertical profiles which are 
not similar. 


The numerical model calculations were made with approximately 20 000 
nodes. Although these are In qualitative agreement with the data, they show 
temperature gradients that are too steep, especially In the transverse direc- 
tion. Underprediction of the mixing was also seen In the single-jet calcula- 
tions of reference 13 where It Is shown that the k-e type of turbulence model 
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underestimates the Intensity. The result In figure 6 Is typical of the numer- 
ical model calculations shown In this paper. 

For the small-orifice case a coarse-grld calculation using less than 6000 
nodes was also performed. This calculation Is shown on the right side of 
figure 6(a) and Illustrates the significant Influence that grid selection can 
have on the solution obtained, and the smearing of the profiles which can occur 
as a result of numerical diffusion. 


Coupled Spacing and Momentum-Flux Ratio 

Examination of the experimental data revealed that similar jet penetration 
Is obtained. Independent of orifice diameter. If orifice spacing and momentum- 
flux ratio are coupled (refs. 2, 5, 12, and 22). For example, low momentum- 
flux ratios require large, widely spaced holes, whereas smaller, closely spaced 
holes are appropriate for high-momentum flux ratios, as shown In figure 7. The 
duct cross section Is shown to the right of the three-dimensional oblique and 
Isotherm contour plots for each configuration. It follows that for low 
momentum-flux ratios (large spacing) a greater axial distance Is required for 
equivalent mixing. 

In general, jet penetration and centerplane profiles are similar when the 
orifice spacing and the square root of the momentum-flux ratio are Inversely 
proportional; that Is, 


C = (S/H q ) y/J (3) 

For single-side Injection, the centerplane profiles are approximately 
centered across the duct height and approach an Isothermal distribution In the 
minimum downstream distance when C = 2.5. This appears to be Independent of 
orifice diameter, as shown In both the calculated and experimental profiles In 
figure 8. In equation (3), values of C that are a factor of 2 or more smal- 
ler or larger than the optimum values correspond to underpenetration or over- 
penetration, respectively. (Figs. 5 and 6 and table I). A summary of the 
spacing and momentum-flux ratio relationships for single-side Injection Is 
given In table II. 


Variable Temperature Main Stream 

The Influence of a nonlsothermal main-stream flow on measured profiles 
for Intermediate momentum-flux ratios with S/Hq » 0.5 and Hq/D * 4 Is 
shown In figure 9. The corresponding Isothermal main-stream case Is shown In 
the top row. In the center row of the figure, the upstream profile (left 
frame) Is coldest near the Injection wall, whereas In the bottom row, the 
upstream profile (left frame) Is coldest near the opposite wall. For the 
definition of e In this figure, T m Is the hottest temperature In the 
main stream for each case. 

Experimental, empirical, and numerical results for the top-cold case are 
shown In figure 10. The empirical calculations are from a superposition of the 
upstream profile and the corresponding jets-ln-an-lsotherfnal-maln-stream dis- 
tribution (ref. 22). Although this gives a good first-order approximation. It 
should be noted that with a variable temperature main stream there can be 
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cross-stream thermal transport due to the flow of main-stream fluid around the 
jets (and hence, to different y locations), and this Is not accounted for 
when the distributions are superimposed. This becomes apparent If the local 
main-stream temperature T m (y), Is used In the definition of e In 
equation (1). 

In the variable temperature main- stream case, the numerical model results 
agree well with the experimental data, especially on the jet centerplane, but 
the transverse mixing Is underpredicted, as In the corresponding Isothermal 
main-stream case shown In figure 8(b). 


Slanted Slots 

Figure 11 shows experimental and calculated three-dimensional oblique 
e distributions for slanted slots at Intermediate momentum-flux ratios. These 
slots had an aspect ratio (length/width) of 2.8, with their major axes slanted 
at 45° to the main-stream flow direction. The orifice spacing and size are 
S/Hg = 0.5 and Hg/D = 4, respectively. The temperature distributions In this 
figure may be compared with those for equivalent-area circular holes In fig- 
ure 8(b). As noted In reference 23, the penetration and mixing of jets from 
the slanted slots are noticeably less than those of jets from the circular 
holes. The normally symmetric vortex pair Is asymmetric In this case, as Is 
apparent In the experimental profiles In figure 11 (and In the figures In 
ref. 23). These profiles also show that the centerplanes of the jets shift 
laterally with Increasing downstream distance. 

The empirical model calculations Include a modification to account for the 
observed centerplane shift, but they do not model the asymmetry (refs. 19 
and 23). The numerical calculations for this case exhibit both the centerplane 
shift and the asymmetry, and they are good In the context of the qualitative 
agreement seen throughout the comparisons given In this paper. 


Double Rows of Hoi es 

Figure 12 shows experimental and calculated temperature distributions for 
an orifice plate with two Inline rows of jets (S x /Hg = 0.5) from circular 
orifices. It was observed from the experimental profiles In reference 23 that 
the two configurations have very similar temperature distributions, and this 
Is seen In the calculated profiles as well. In this case the empirical model 
calculations are derived by superimposing the distributions from the two rows. 

Both experimental and calculated temperature distributions are shown In 
figure 13 for a double-row configuration when S x /Flg = 0.25 and the trailing 
row has twice as many orifices as the lead row. Note that the orifice area Is 
the same for both rows. These profiles show the dominance of the lead row In 
establishing the jet penetration and first-order profile shape (ref. 23). As 
with the double row of Inline holes, the empirical calculations for this case 
were obtained by superimposing separate calculations for the two rows. 
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Opposing Rows of Jets 


The remainder of this paper presents results for two-side Injection from 
opposing rows of Jets with (1) the top and bottom Jet centerlines directly 
opposite each other and (2) the top and bottom Jet centerlines staggered In the 
z direction. The experimental results are shown and compared with the single- 
side results In figures 14 and 16. In these figures, a duct cross section, 
drawn to scale. Is to the left of the data. 


Opposed Rows of Inline Jets 

Figure 14 shows a comparison at Intermediate momentum-flux ratios for 
single-side and opposed-jet Injection. For these momentum-flux ratios, an 
appropriate orlf Ice-spaclng-to-duct-helght ratio for optimum single-side mixing 
Is approximately 0.5 (see eq.(2)), as confirmed by the profiles In figure 7. 

For opposed-jet Injection, with equal momentum-flux ratios on both sides, 
the effective mixing height Is half the duct height because, as reference 3 
shows, the effect of an opposite wall Is similar to that of the plane of sym- 
metry In an opposed-jet configuration. Thus, the appropriate orlf Ice-spaclng- 
to-duct-helght ratio for opposed-jet Injection at these Intermediate 
momentum-flux ratios would be about S/Hq = 0.25. Dimensionless temperature 
distributions downstream of Jets with this spacing are shown In the bottom row 
of figure 14, and the two streams do Indeed mix very rapidly. Note that since 
the orifices In figures 14(a) and (b) are the same size, the Jet to main-stream 
flow ratio Is four times greater for opposed-jet Injection than for single-side 
Injection. If It Is desired to maintain an equal flow rate, the orifice diam- 
eter must be halved, since there Is Injection from both sides and opposed-jet 
Injection requires twice as many holes In the row as optimum single-side 
Injection. 

Figure 15 shows experimental and calculated profiles for opposed rows of 
jets with Identical orifice spacing and diameter and with the orifice center- 
lines In line. The empirical model predicts opposed-jet Injection very well, 
as the experimental profiles on both sides of the plane of symmetry support the 
Gaussian profile assumption. The penetration and profile shape calculated with 
the numerical model are In good agreement with the data, but the mixing Is 
otherwise underpredicted, as evidenced by the steep transverse and lateral 
gradients seen In almost all the previous calculations also. 


Opposed Rows of Staggered Jets 

Figure 16 shows comparisons of single-side and staggered jet Injection 
for Intermediate momentum-flux ratios. Since the effective mixing height for 
opposed Inline Injection was half the duct height. It was assumed that the 
effective orifice spacing for staggered jets would be half the actual spacing. 
Thus, to maintain an optimum coupling of the effective spacing and the 
momentum-flux ratio, the orifice spacing for opposed staggered configurations 
should be double that which Is appropriate for single-side Injection. 

This hypothesis Is verified by the rapid mixing of the two streams In the 
bottom row of profiles In figure 16. This figure shows clearly that a conflg- 
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uratlon which mixes well with one-side Injection performs even better when 
every other orifice Is moved to the opposite wall. 

Empirical and numerical model calculations for an opposed row of staggered 
jets are compared with the data In figure 17. The empirical model does not 
handle this complex case well, as the fluid dynamic Interactions here are not 
amenable to a direct extension of the simple Gaussian profile and superposition 
type of modeling appropriate for most of the single-side and opposed-jet cases 
of Interest. The numerical model calculations are not In appreciably better 
agreement with the data than are the empirical model results, however, because 
the mixing Is underpredicted here as In the previous cases. 

A summary of the spacing and momentum-flux ratio relationships which give 
optimum mixing for opposed rows of Inline and staggered jets Is given In 
table II. 


LIMITATIONS AND APPLICABILITY 
Empl rlcal 

The empirical model results presented here show that correlating experi- 
mental data can provide an excellent predictive capability within the parameter 
range of the generating experiments, provided the experimental results are 
consistent with the assumptions made In the empirical model. These models 
must, however, be used with caution, or not at all, outside this range. 

The ranges of the experimental variables on which this empirical model was 
based are given In table III. The density ratio, momentum-flux ratio, orifice 
spacing, and orifice size were the primary Independent variables. This table 
also gives ratios that are derived from the primary variables: the orifice to 

main-stream area ratio, the jet-to-total mass flow split, and the parameter 
coupling the spacing and momentum-flux ratio. Not all combinations of the 
primary variables In the table were tested; only those combinations which are 
within the range given for the derived variables represent conditions that are 
within the range of the experiments. 

Examining the results In figures 5 to 17 In the context of equation (3) 
suggests that, In general, the empirical model provides good temperature field 
predictions for single-side Injection when 1 < C < 5. Similarly, good predic- 
tions are obtained for opposed Inline jets provided that 0.5 < C < 2.5. This 
model does not work well for Impinging flows because the experimental tempera- 
ture distributions are not consistent with the assumption of Gaussian profile 
similarity In the empirical model. The experimental profiles for conditions 
giving optimum mixing In opposed staggered-jet configurations are also somewhat 
at variance with the model assumptions; In these cases, satisfactory agreement 
with the data must be considered fortuitous. 

A major weakness of the empirical model used here In (refs. 17 to 19) and 
In previous versions (refs. 4 and 5) Is that the form of the correlations pre- 
cludes their use for semi-confined flows (large Hg/D or S/D), single-jet 
flows, or flows In which It Is known a priori that the primary assumptions In 
the model will be Invalid. 
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Numerical 


The numerical model Is not subject to the Inherent limitation of the 
empirical model regarding profile shape and confinement. Thus, three- 
dimensional codes can provide calculations for complex flows for which the 
assumptions In the empirical model are known to be Invalid or outside the range 
of available experiments. Furthermore, numerical models provide calculations 
for all flow field parameters of Interest, not just those that happen to have 
been empirically correlated. 

The numerical calculations correctly show the trends which result from 
variation of the Independent flow and geometric variables, although the results 
consistently exhibit too little mixing. The numerical model calculations for 
the slanted slots and staggered jet cases are encouraging because the experi- 
mental data for these cases show profiles that are not consistent with the 
primary assumptions In the empirical model. 

The numerical calculations performed are shown to be grid sensitive, and 
false diffusion Is known to be present. Uncertainties also exist In these 
calculations regarding the validity of turbulence model assumptions and due to 
unmeasured (and hence assumed) boundary conditions. The results shown here are 
not Intended to represent the best agreement possible from numerical models at 
this time. Better temperature field agreement could undoubtedly have been 
achieved by adjusting model constants and/or Inlet boundary conditions. But, 
since this was not necessary to satisfy the present objective of evaluating the 
potential of these codes vis-a-vis combustor dilution zone flow fields, and 
because the mean temperature was the only parameter compared, no adjustments 
were made. 

Thus, as with previous assessments In references 14 to 16, three- 
dimensional calculations, such as those In this paper, should be considered as 
only qualitatively accurate at this time, and three-dimensional codes of this 
type are useful primarily In guiding design changes or In perturbation analy- 
ses. The three-dimensional code used herein, although sufficiently promising 
to justify further development and assessment. Is not a practical tool for 
general engineering use In Its present form. Codes with Improved numerics, 
accuracy, and turbulence models should provide more quantitative predictions. 


CONCLUSIONS 

The principal conclusions from the experimental results reviewed herein 
are as follows: 

1. Variations of momentum-flux ratio and of orifice size and spacing have 
a significant effect on the flow distribution. 

2 . Similar distributions can be obtained. Independent of orifice diameter, 
when momentum-flux ratio and orifice spacing are coupled. 

3. A first-order approximation of the mixing of jets with a variable- 
temperature main stream can be achieved by superimposing the jets-ln-an- 
Isothermal-maln-stream and upstream profiles. 


159 



4. The penetration and mixing of jets Issuing from 45° slanted slots are 
less than those from equivalent-area circular holes. 

5. With the same orifice spacing In (at least) the lead row, double rows 
of jets have temperature distributions similar to those from a single row of 
equally spaced, equivalent-area circular orifices. 

6. For opposed rows of jets, with the orifice centerlines In line, the 
optimum ratio of orifice spacing to duct height Is one-half the optimum value 
for single-side Injection at the same momentum-flux ratio. 

7. For opposed rows of jets with the orifice centerlines staggered, the 
optimum ratio of orifice spacing to duct height Is twice the optimum value for 
single-side Injection at the same momentum-flux ratio. 

Temperature field measurements from the experiments cited previously are 
compared with distributions calculated using an empirical model based on 
assumed vertical profile similarity and superposition and with calculations 
made using a three-dimensional elliptic code with a standard k-e turbulence 
model. The results can be summarized as follows: 

Empirical model calculations provide very good results for modeled param- 
eters within the range of experiments whenever the primary assumptions In the 
model are satisfied. 

Three-dimensional code calculations made In this study correctly approxi- 
mate the trends which result from varying the Independent flow and geometric 
variables, but they consistently exhibit too little mixing. (The advantage of 
these models Is that they can predict all flow field quantities, flows outside 
the range of experiments, or flows where empirical assumptions are Invalid.) 
Numerical calculations should yield more quantitative predictions with 
Improvements In numerics, accuracy, and turbulence models. 
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TABLE I. - FLOW AND GE0ME1RY CONDITIONS 


Figure 

S/H 0 

V“ 

A. / A 
J m 

C d 

DR 

Vi 

V“t 

c a 

2 

0.5 

4 

0.10 

0.76 

2.2 

26.2 

0.36 

2.56 

5(a), 6(a) 

.25 

8 

.05 

.60 

2.1 

22.4 

.17 

1.18 

5(b), 6(b) 

1 .0 

4 

.05 

.67 

2.2 

23.5 

.19 

4.85 

7(a) 

1 .0 

4 

.05 

.73 

2.1 

5.3 

.11 

2.30 

7(b), 14(a) 

.5 

8 

.025 

.61 

2.2 

28.4 

.11 

2.66 

7(c) 

.25 

8 

.05 

.61 

2.3 

92.7 

.30 

2.60 

8(a) 

.5 

5.7 

.05 

.71 

2.2 

25.4 

.21 

2.52 

8(b), 9(a), 16(a) 

.5 

4 

.10 

.61 

2.1 

18.6 

.27 

2.16 

t>9( b) , 10 

.5 

4 

.10 

.61 

1 .8 

31 .3 

.31 

2.80 

C9(C) 

.5 

4 

.10 

.68 

2.2 

24.4 

.31 

2.47 

dll 

.5 

4 

.10 

.66 i 

2.2 

27.1 

.33 

2.60 

12 

.5 

5.7 

.05 

.65 

2.2 

26.3 

.33 

2.56 


.5 

5.7 

.05 

.66 

2.2 

26.9 

— 

2.59 

13 

.5 

5.7 

.05 

.69 

2.2 

26.8 

.34 

2.59 


.25 

8 

.05 

.70 

2.2 

26.6 

- 

1 .29 

14(b), 15 

.25 

8 

.10 

.65 

2.1 

25.0 

.32 

1 .25 

16(b), 17 

1 .0 

4 

.10 

.65 

2.1 . 

27.6 

.33 

5.25 


ar = ^s/HaW/T 

. - - "O' V ' - 

t>Top cold. 
c Top hot. 

d45° slanted slots. 
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TABLE II. - SPACING AND MOMENTUM-FLUX RATIO 
RELATIONS 


Configuration 

C = (S/Ho) V^T 

Single-side Injection 


Underpenetration 

<1.25 

Optimum 

2.5 

Overpenetration 

>5 

Opposed rows of jets 


Inline optimum 

1 .25 

Staggered optimum 

5 


TABLE III. - RANGE OF INDEPENDENT FLOW 
AND GEOMETRIC VARIABLES INVESTI- 
GATED IN REFERENCES 17 TO 19 


DR 

... 0.5 to 2.5 

J 

.... 5 to 105 

S/Hq 

. . . 0.125 to 1 

Ho/D 

4 to 16 

Aj/A m 

. . 0.025 to 0.1 

Wj/WT ....... 

. . 0.075 to 0.33 

C = (S/Ho)^T . . 

. . . . 0.5 to 10 
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Figure 3. - Dilution-jet mixing orifice configurations. 



Figure 4. - Schematic of typical vertical 
temperature profile showing scaling 
parameters in empirical model. 
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(b) S/H 0 1 1; Hg/D ■ 4; J *23.5. 

Figure 6. - Effect of varying orifice spacing at constant area on measured and calculated temperature distributions when x/Hq * 0.5 (AJA n 
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EXPERIMENTAL EMPIRICAL NUMERICAL 

35 x 33 x 17 



(a) Hq/D =5.7; Aj/A m =0.05; J =25.5. 


45 x 23 x 19 



0 1 0 1 0 1 

e 6 e 

(b) Hq/D = 4; Aj/A m = 0. 10; J = 18. 6. 


Figure 8. - Effect of varying orifice diameter at constant spacing on measured and calculated 
temperature distributions when x/H 0 = 0. 5 (S/Hq = 0. 5). 


x/H 0 = 0. 5 1 2 



(a) Isothermal mainstream; J = 18.6. 



(c) Top hot; J = 24.4. 


Figure 9. - Influence of nonisothermal mainstream on measured temperature profiles (S/Hq = 0.5; Hq/D = 4, Aj/A m c 0. 10). 
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EXPERIMENTAL EMPIRICAL NUMERICAL 

45 x 23 x 19 
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(b) x/H 0 * 1. 0. 

Figure 10. - Measured and calculated temperature distribution for jets injected into a nonisothermal 
mainstream; top cold (S/Hg * 0. 5, Hq/D * 4, J ■ 31. 3). 




0 0 0 


(c) x/H 0 ■ 1. 0 

Figure 11. - Measured and calculated temperature distributions for slanted slots at an intermediate momentum-flux ratio (S/H 0 ■ 0. 5, 
Hq/D ■ 4, J ■ 27. 1). 
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e e e 

(b) Opposed row (inline) injection: S/Hg ■ 0.25; Aj/A m - 0. 10; J * 25. 0. 

Figure 14. - Comparison between single-side and opposed- jet injection (Hg/D - 8). 


EXPERIMENTAL EMPIRICAL NUMERICAL 


45 x 26 x 17 



0 1 0 1 0 1 


e e e 

(b) x/H 0 = 1. o. 

Figure 15. - Measured and calculated temperature distributions for opposed rows of inline jets (S/Hg ■ 0.25, Hg/D - 8. 
Aj/A m * 0. 10, J « 25). 
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(a) Single-side (top) injection: S/Hq = 0.5; J ■ 18.6. 



0 0 0 
(b) Opposed row (staggered) injection: S/Hq = 1 ; J = 27.6. 

Figure 16. - Comparison between single-side and staggered jet injection (Hq/D = 4, Aj7A m « 0. 10). 


EXPERIMENTAL EMPIRICAL NUMERICAL 

22 x 27 x 33 



0 0 6 


(b) x/H q = 1. 0. 

Figure 17. - Measured and calculated temperature distributions for opposed rows of staggered jets (S/H 0 = 1. Hq/D = 4, 
Aj/A m =0.10, J =27.6). 
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THEORETICAL KINETIC COMPUTATIONS IN COMPLEX REACTING SYSTEMS 

David A. Blttker 
NASA Lewis Research Center 
Cleveland, Ohio 


This paper describes NASA Lewis' studies of complex reacting systems at 
high temperature. The changes which occur are the result of many different 
chemical reactions occurring at the same time. Both an experimental and a 
theoretical approach are needed to fully understand what happens In these sys- 
tems. The latter approach Is discussed herein. We present the differential 
equations which describe the chemical and thermodynamic changes, and we 
describe their solution by numerical techniques using a detailed chemical 
mechanism. Several different comparisons of computed results with experimen- 
tal measurements are also given. These Include the computation of (1) species 
concentration profiles In batch and flow reactions, (2) rocket performance In 
nozzle expansions, and (3) pressure versus time profiles In hydrocarbon Igni- 
tion processes. The examples Illustrate the use of detailed kinetic computa- 
tions to elucidate a chemical mechanism and to compute practical quantities 
such as rocket performance. Ignition delay times, and Ignition lengths In flow 
processes. 


INTRODUCTION 

For many years Lewis has been studying complex, reacting gas-phase sys- 
tems at high temperature. The changes that occur In such systems are the 
result of many Individual chemical reactions occurring at the same time. We 
use a two-step approach to understand what happens In such a system: (1) 

experimental measurement of the temporal changes In temperature, pressure, or 
composition of the system, and (2) theoretical computation of these changes In 
an attempt to match the experimental results. This report describes our work 
In the computation of complex-system chemical kinetics as focused on combus- 
tion systems. Brabbs describes our experimental efforts In his paper for this 
symposium. Early attempts to compute the progress of complex reactions 
Involved various simplifications. These were of two types: (1) one or two 
global reactions were substituted for the actual set of many Individual, 
simultaneous reactions which occur In a gaseous system with several species 
present, and (2) reactions Involving some of the very reactive atoms and radi- 
cals were assumed to be very fast and, therefore, always In chemical equilib- 
rium. Differential equations for the rates of change of a few key species 
were solved analytically. Although both of these approaches have had limited 
success, they only explain some of the general features of a complex reaction. 
They are, at best, simplified approximations of the actual process, which are 
only valid under a limited set of conditions. The most fruitful approach to 
complex kinetics computations Is to numerically Integrate a system of differ- 
ential equations derived from the laws of conservation of mass, energy, and 
(If flow Is Involved) momentum. These differential equations Involve the 
rates of change of species concentrations oj (moles per unit mass of mix- 
ture) and temperature T, and they may also Involve density p and velocity 
V for a flow process. When these differential equations are solved, a set of 
Individual reactions Is assumed to occur simultaneously among all the species, 
and these reactions are usually assumed to be reversible (with the ratio of 
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the forward to reverse rate constant being equal to the equilibrium constant 
for any reaction) . 

There are two ways to use these detailed kinetic computations. First, 
they may be used to elucidate a reaction mechanism or to determine a single 
unknown rate constant In conjunction with experimental data. On the other 
hand, they may be used with a known mechanism to compute practical combustion 
quantities such as 

(1) Kinetic rocket or jet engine performance 

(2) Composition of engine exhaust gases 

(3) Ignition delays In combustible mixtures 

(4) Ignition lengths and nozzle performance in supersonic combustion 

Examples of some of these computations performed at NASA Lewis are given later 
In this paper. 


GENERAL THEORY 

A reacting system containing NS species Is assumed to have NR Inde- 
pendent chemical reactions proceeding simultaneously. The general set of 
reactions Is written 

NS 


1=1 

Using the law of mass action, we can write the forward and reverse rate of 
each reaction as 


NS 



where p Is the mass density of the mixture and the concentration of 
species Si, (moles per unit mass of mixture). The rate constant for the for- 
ward reaction kj Is written In the modified Arrhenius form 

n -E./RT 

kj = AjT J e 3 (3) 




j=l ,2, . . . , NR 
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The reverse rate constant k_j 
reverslbl lity: 


In these last equations, 

Ej = activation energy for the 
T = temperature, K 
R = universal gas constant 
(K e q)j = equilibrium constant for reaction j 

For any reacting system, batch or flow, the appropriate continuity equa- 
tions are differentiated to give a set of differential equations for the vari- 
ation of the concentrations and temperature with time. For a flow reaction, 
differential equations for density and velocity will also be obtained. The 
resulting system of equations Is solved numerically In combination with the 
Ideal gas equation of state: 


Is obtained from the law of microscopic 


k _j * K eq*j 


( 4 ) 


jth reaction, energy/mol 


P 

p = NS 

XT £ «, 
1=1 


( 5 ) 


where p Is the pressure of the mixture. 

For the case of constant pressure and adiabatic batch reaction, the equa- 
tions are as follows: 


d „ w 

dt - "7 U1 ' 2 " S 

where W-j , the molar rate of formation of species , Is given by 

NR 

W 1 ■ L ( *ij ~ ”13 )(R 3 " R -3’ 

3=1 

and 


dl 

dt “ ” NS 


NS 

w i h i 


E cr. C , 

i-i 1 Pr1 


( 6 ) 


( 7 ) 


( 8 ) 
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where 


= molar enthalpy of species 
c p> i = molar heat capacity of species 

The standard Runge-Kutta or other explicit numerical techniques cannot be used 
to solve these equations. Although the differential equations are Inherently 
stable, they have widely different time constants for relaxation to their 
final equilibrium solution. Therefore, the step sizes required for solution 
stability are prohibitively small. These systems of equations are called 
stiff, and the problem of solving them accurately has been studied for many 
years. Several new Implicit Integration methods have been developed for solv- 
ing these stiff equations. Although they have been quite successful In 
Improving the accuracy and efficiency of the solutions, there Is still a need 
for additional development of better methods. Problems of accuracy and of 
efficient step-size selection can still arise. Methods need to be developed 
for automatic selection of the optimum Integration parameters. Therefore, 
research Is actively proceeding at many Institutions to develop more efficient 
Integration methods (refs. 1 to 4). 

There Is another kinetic combustion model of Interest to practical engine 
designers - namely, highly backmlxed reacting flow In what Is usually called 
the well-stirred reactor. We assume the limiting condition of zero-dimensional 
flow or Instantaneous backmlxlng of the reacted gases with the cold, unreacted 
gases. Although this model Is an oversimplification of highly turbulent 
reacting flow, It Is a very useful first-order approximation for some practi- 
cal reactlng-flow systems. 

The process Is a constant pressure combustion with mass flow rate m 
through a reactor of constant volume v. The average residence time In the 
reactor Is 


For each species we can write a continuity relation which says the fol- 
lowing: The difference between the species molar flow rate Into and out of 

the reactor Is equal to Its rate of formation (or destruction) by chemical 
reaction In the reactor. These continuity equations are 

NR 

® <°1 - «i) + XI (v ij " v 1j )(R j ' R -j } 1=1 ’ 2 NS (10) 

The following energy conservation equation can also be written (If we assume 
that the process Is adiabatic): 

NS 

E (o^ - o*h*) =0 (11) 

In these equations * Indicates the unreacted gas mixture. Equations 
(10) and (11) are a set of nonlinear, algebraic equations for and the 
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reactor temperature. They are solved by the Newton-Raphson Iterative proce- 
dure. Logarithmic Increments of the variables are used to avoid numerical 
problems (ref. 5). 


NASA GENERAL CHEMICAL-KINETICS CODE 

Many computer codes have been published which perform different kinds of 
chemical-kinetic computations. At NASA Lewis we published a general chemical- 
kinetics code, GCKP84 (ref. 6). It performs a wide variety of chemical- 
kinetics computations with convenience and efficiency. 

It Is designed to perform chemical-kinetics computations for several 
reaction models. Including the following: 

(1) General reaction In either a batch system or a one-dimensional, 
frlctlonless plug flow 

(2) Combustion reaction In a well-stirred reactor (highly backmlxed flow) 

(3) Reaction behind a shockwave with boundary layer corrections 

(4) Ignition processes In either a batch or flow system 

(5) Nozzle expansion reactions 

For each of these models the following general features of the code apply: 

(1) Any chemical system may be used for which reaction rate constant data 
and species thermodynamic data are known. 

(2) The process may be adiabatic, or the heat transfer between the 
reaction and Its environment may be considered. 

(3) A new efficient Integration technique Is employed (ref. 7). 

(4) Any chemical reaction of the form aA + bB ^dD + eE may be used, 
Including photochemical and Ionic reactions. 

(5) Simplified Input for combustion reactions may be used. 

The code contains an option to compute rocket performance parameters for 
nozzle flow and can be conveniently used to compute Ignition lengths In super- 
sonic flow. An option Is also provided to perform a well-stirred reactor 
computation and, then. Immediately expand the products through an assigned 
nozzle profile. 


EXAMPLES OF CHEMICAL KINETIC COMPUTATIONS WITH GCKP84 
Hydrogen-Oxygen Batch Reaction at Constant Temperature 

Measurements of hydrogen peroxide formation In the constant-volume. Iso- 
thermal reaction of hydrogen and oxygen at 500 °C were reported many years ago 
by Baldwin (ref. 8). The Important reactions In the hydrogen-oxygen mechanism 
are as follows: 


H 2 + 0 2 = OH + OH 
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H 2 + OH = H 2 0 + H 
H + 0 2 = OH + 0 
0 + H 2 = OH + H 
*H + 0 2 + M = H0 2 + M 
H + H0 2 = OH + OH 
*H 2 0 2 + M = OH + OH + M 
*ho 2 + H0 2 = H 2 0 2 + o 2 
*H 2 + H0 2 = H 2 0 2 + H 

The rate constants for many of these reactions have been measured fairly 
accurately, but some of the reactions Involving H 2 0 2 and H0 2 still have a sig- 
nificant uncertainty. The reactions Important In determining H 2 0 2 concentra- 
tion are Indicated by an asterisk. GCKP84 was used to compute the H 2 0 2 versus 
time profile. After making several variations of the rate constants for the 
two most uncertain and Important reactions (H 2 + H0 2 H 2 0 2 + H and H0 2 + 

H0 2 ^ H 2 0 2 + 0 2 ), we obtained the agreement between the experimental and com- 
puted results shown In figure 1. This example Illustrates the use of kinetics 
computations In obtaining better rate-constant values for one or two key reac- 
tions when all other reactions In a mechanism are fairly well known. Of 
course, we cannot say that we have uniquely determined both reaction rate con- 
stants since we have two adjustable parameters. But we have put limits on 
their uncertainty. Of course, this type of exercise can only be done when the 
experimental variable being matched by computation Is sensitive to small vari- 
ations In the uncertain rate constants. This brings up the following general 
question: For a given change In any rate constant of the reaction mechanism, 

how much does each computed variable change? This Is the relatively new 
field of study In chemical kinetics called sensitivity analysis. 

In practice, even through a chemical mechanism may contain more than 100 
reactions, only the rate constants of 15 or 20 may significantly affect the 
computed results when their values are changed. Recently, methods have been 
developed to systematically compute these effects In the form of sensitivity 
coefficients (refs. 9 to 11). This computation Is performed along with the 
chemical-kinetics computation. The results of the latter computation are used 
as Input to solve a second set of differential equations whose unknowns are 
sensitivity coefficients of the form 


T 



k 0 k 
8. m 


Hm 


i c i 

ak„ 


2 

a c. 


ak ak 

Si m 


( 12 ) 

(13) 


The first-order coefficient effectively gives the percentage change In 

c^ ( 1 th concentration variable) for a given percentage change In the rate 
constant of reaction 8. The second-order coefficient Tia.m gives the effect 
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on of changes In both k a and k m . It is now recognized that a sensi- 

tivity analysis must be performed for any complex reaction In order to pinpoint 
the Important reaction paths and obtain a good understanding of the reaction 
mechanism. 


Computation of Rocket Performance 

One of the first uses of chemical kinetic computations was In the compu- 
tation of rocket performance using a one-dimensional flow model. This per- 
formance Is measured by the velocity of the hot gas at the nozzle exit point. 

As the gas expands and cools. Its velocity (i.e., kinetic energy) Is kept as 
high as possible by chemical recombination reactions which occur In the noz- 
zle. These reactions convert high potential-energy atoms and radicals Into 
low-energy stable molecules. Maximum performance Is obtained If these reac- 
tions maintain equilibrium conditions in the very short nozzle-residence time. 
But this Is often not the case since the finite rates of the reactions are not 
fast enough. By knowing the rate constants for the Important recombination 
processes one can compute kinetic performance that more accurately and realis- 
tically reflects the experimentally measured performance. An example of this 
is shown In figure 2 for the oxygen difluoride - dlborane system (ref. 12). 

The three curves of performance as a function of oxidant/fuel ratio (fig. 2(b)) 
show the maximum and minimum theoretical performance (assuming equilibrium and 
frozen conditions) as well as the klnetlcally limited performance. The latter 
agrees rather well with actual measurements corrected for various losses. The 
reaction mechanism used Is shown in table I. This system Is Interesting 
because only one reaction rate constant has to be accurately known In order to 
compute the kinetic curve. This is the recombination of hydrogen atoms to 
give molecular hydrogen. This fact was determined by a simple sensitivity 
analysis. Computed performance is unchanged for wide variations of the other 
rate constants. The reason for this situation can be seen from the heat 
release data In table I. The hydrogen atom recombination accounts for well 
over 50 percent of the heat release in the process, and this, of course, con- 
trols performance. 


Two-Stage Well-Stirred Reactor 

The practical use of the well- stirred reactor model of GCKP84 is Illus- 
trated by Its simulation of a two-stage turbulent combustor. A theoretical 
and experimental study of two-stage, rich- lean combustion as a method of con- 
trolling nitrogen oxide emissions was performed at Lewis. This technique was 
suggested as a way of reducing the conversion of the organic nitrogen compounds 
in hydrocarbons (fuel-bound nitrogen) to nitric oxide during combustion. The 
objective of the work was to study the effect of operating conditions on emis- 
sions of both nitrogen oxides (NO x ) and carbon monoxide (CO) (ref. 13). 

Both the nitrogen and the hydrogen content of the fuel were changed to simu- 
late coal-derived syncrude fuels. These fuels contain high percentages of 
aromatic hydrocarbons and therefore have a lower hydrogen-to-carbon ratio than 
today's petroleum-based fuels. Experiments were performed in the two-stage 
flame-tube apparatus shown In figure 3. The fuels were mixtures of propane, 
toluene, and pyridine blended to give a range of hydrogen and nitrogen con- 
tents. The primary-zone equivalence ratio was varied from 0.7 to about 1.8. 

The products of first-stage combustion were rapidly diluted with air Injection 
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and burned at a final equivalence ratio of 0.5. The final N0 X and CO concen- 
trations were measured. The two-stage flame tube was modeled analytically by 
a two-stage stlrred-reactor computation closely simulating the experimental 
fuels. Instantaneous mixing of the secondary dilution air was assumed. A 
detailed chemical mechanism was used for propane and toluene oxidation and for 
N0 X formation. The fuel-bound nitrogen was Input as nitrogen atoms and all 
mixing and heat-transfer effects were neglected. This very simplified model 
of a quite complex process was able to predict most of the observed trends of 
N0 X and CO formation with variation of nitrogen and hydrogen concentrations 
In the fuel. The model's prediction of other trends, which were not measured 
experimentally, gave additional Information about the advantages and disadvan- 
tages of using this rich-lean combustion for emissions control. A comparison 
of some experimental and computed results for N0 X and CO concentration as a 
function of primary equivalence ratio Is given In figures 4 and 5. The simple 
model qualitatively predicts the observed trends. It can be seen that using a 
rich primary equivalence ratio does reduce N0 X formation very significantly. 
However, at the same time, this technique Increases the formation of CO, which 
Is also undesirable. Thus, one set of conditions cannot minimize both pollut- 
ants, and tradeoffs will have to be made to obtain desired emissions levels 
for any practical situation. This work has shown that a simple kinetic model 
can sometimes be used to qualitatively explain the observed trends In a highly 
complicated combustion system. 


Ignition Delays In Benzene-Oxygen-Argon Mixtures 

Our current work Involves the study of the mechanism of hydrocarbon oxida- 
tion. We are studying the aliphatic hydrocarbons by measuring Ignition delays 
and concentration profiles during the shock-heated oxidation of hydrocarbon- 
oxygen-argon mixtures. Then a detailed chemical mechanism Is formulated and 
used In detailed kinetics computations to match the observed experimental data 
over a wide range of Initial composition, temperature, and pressure. We are 
also studying the simplest aromatic, benzene. A partial mechanism for benzene 
oxidation Is as follows: 


c 6 h 6 + °2 = c 6 h 5 + h0 2 
C 6 H 6 = C 6 H 5 + H 

c 6 h 6 = c 4 h 4 + c 2 h 2 
H + 0 2 = OH + 0 

H + 0 2 + M = H0 2 ♦ M 
0 + C 6 H 6 = C 6 H 5 + OH 
OH + C 6 H 6 = C 6 H 5 + H 2 0 
C 6 H 5 + ° 2 = C 6 H 5 0 ♦ 0 
C 6 H 5 0 + M = C 5 H 5 + CO + M 
CO + OH = C0 2 + H 
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The experimental program Is described by T.A. Brabb's paper presented at this 
symposium. I attempted to match experimental Ignition delay data for benzene- 
oxygen-argon mixtures that were taken by another Investigator at Lewis. The 
reaction zone behind the experimental reflected shock Is approximated by a 
constant volume batch reaction. The theoretically computed temperature and 
pressure (for nonreacting conditions) behind the shock are used as starting 
conditions. A typical, computed profile of pressure as a function of time Is 
shown In figure 6. The Ignition delay t Is determined from the first sig- 
nificant pressure rise; this Is similar to determining the experimental t 
value from the measured pressure trace. The mechanism we used Includes over 
100 reactions Involving benzene and Its degradation products, acetylene, phe- 
nol, phenylacetylene, ethylene, and methane. It Includes all the reactions of 
the hydrogen-oxygen system as well. Values from the literature were used for 
all reaction rate constants except the quite uncertain ones of the benzene, 
phenol, and phenylacetylene reactions. These were varied within reasonable 
limits to get the best overall agreement between computed and experimental 
t values. Starting mixture equivalence ratios ranged from 0.5 to 2, tempera- 
tures ranged from about 1200 to 1700 K, and Initial pressures ranged from 2 to 
6 atm. A simple sensitivity study made by changing rate constants one at a 
time has shown that the x values are most sensitive to the rate constant for 
the C5H5 + O2 reaction, which Is quite uncertain. Other rate constants which 
are uncertain and have a significant effect are those for C6H6 * C4H4 + C2H2 
and C&H5 C4H3 + C2H2. In figure 7 we show some comparisons between com- 
puted and experimental Ignition delays. We have used only the experimental 
data points which are considered the most accurate for the mechanism matching. 
These are the data with x values of 100 msec or greater. Shorter Ignition 
delays may have significant error due to the nonuniformity of the reaction 
mixture and the heating effects. The computed line for log x versus 1/T 
comes within ±50 percent for all but one point which differs by about a factor 
of two. While this mechanism does fairly well. It Is certainly not complete 
and has two significant problems. First, It tends to predict delays longer 
than the experimental Ignition delays for the lean mixtures and shorter than 
the experimental delays for the rich mixture. Second, It does not predict the 
experimentally observed effect of dilution with argon at a constant equiv- 
alence ratio of 1.0. Therefore, more work needs to be done on this oxidation 
mechanism. A full understanding of the oxidation will probably not be achieved 
until further experimental data are obtained on concentration profiles of a key 
species, such as the phenyl radical or carbon dioxide. Work along these lines 
Is proceeding at Lewis and other laboratories. 
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TABLE I. - 0F 2 -B 2 H 6 MECHANISM 
[Equivalence ratio = 1.63.] 


Reaction 

Net reaction 
rate, 
percent 

Heat of 
reaction, 
AH, 

kcal/mole 

H + H + Mt±H 2 + M 

32.3 

-104 

B0 + HF ^ BOF + H 

26.8 

-33 

h + f + m^hf + m 

5.1 

-136 

H + OH + M ^ H 2 0 + M 

5.2 

-119 

BF + OH ^ BOF + H 

3.9 

-58 

F + H 2 HF + H 

4.5 

-32 

H 2 + OH ^H20 + H 

9.7 

-15 

BF + 0 BO + F 

9.3 

8 

H + BF 2 r^HF + BF 

2.4 

-25 

BO + F + M ^ BOF + M 

.3 

-169 

BF + 0 + M ^ BOF + M 

.2 

-160 

O + H + M^OH + M 

.3 

-102 

Total 

100.0 



Energy release 
rate, 
percent 


100.0 


H 2 O 2 CONCENTRATION 







EXPERIMENTAL 


COMPUTED 


TIME, sec 

Figure 1. - Computed and experimental results for hydrogen -oxygen reaction. 
Temperature, T, 500 °C. 














(a) Rocket engine conditions for OF 2 -BoH 6 propellant Equivalence ratio, 
1.63. 




2.0 Z4 2. 8 3.2 3.6 4.0 

OXIDANT-FUEL RATIO 


(b) Performance of OFo-B^ rocket propellant Chamber 
pressure, 689 kPa 000 psia). 

Figure 2. - Oxygen difluoride-diborane (OF 2 -B 2 H 5 ) rocket 
combustion. 
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EXPERIMENTAL 

DATA 



11% HYDROGEN 
9% HYDROGEN 
SRC COAL 
SYNCRUDE BLEND, 
8.7% HYDROGEN 
0.95% FBN 

COMPUTED RESULTS 


1.0 1.2 1.4 1.6 

PRIMARY EQUIVALENCE RATIO, <p p 

Figure 4. - Effect of hydrogen content on NO x emissions. Fuel-bound 
nitrogen, 1.0 percent; secondary equivalence ratio, 0.5; secondary 
residence time, = 2 msec. 


PERCENT 



PRIMARY EQUIVALENCE RATIO, <p p 

Figure 5. - Effect of percent hydrogen on carbon monoxide 
emissions for hydrocarbon combustion. Fuel-bound ni- 
trogen, 0 percent; secondary equivalence ratio, 0. 5 ; 
secondary residence time, s 2 msec. 
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4.0 

3.6 


3.2 



TIME, Msec 


Figure 6. - Computed pressure as function of time for benzene-oxygen-argon shock ignition. Initial 
temperature, 1366 K; equivalence ratio, <p, 2.0. 


EQUIVALENCE 

RATIO, 

<P 


□ 0.5 

2000 h A 1.0 



Figure 7. - Ignition delay data for benzene-oxygen-argon mixtures. 
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EXPERIMENTAL INVESTIGATION OF PILOTED FLAMEHOLDERS 


C.F. Guo, Y.H. Zhang, and Q.M. Xle 
Gas Turbine Research Institute 
Jlangyou County of Sichuan, People's Republic of China 


Four configurations of piloted flameholders were tested. The range of 
flame stabilization, flame propagation, pressure oscillation during Ignition, 
and pressure drop of the configurations were determined. Some tests showed a 
very strong effect of Inlet flow velocity profile and flameholder geometry on 
flame stabilization. 


INTRODUCTION 

Pressure oscillation during the transient period from Ignition to maximum 
augmentation Is a critical problem for the turbofan augmentor. Since this 
pressure oscillation may cause fan flow to stall or surge. It Is Important to 
know the peak pressure during augmentor Ignition and to take some kind of 
technical measure to reduce this peak pressure. Such a procedure Is called 
"soft Ignition." 

The combustion process In the fan flow of a turbofan augmentor Is under a 
severe condition. It Is helpful to understand how to make use of the hot com- 
bustion products In the core stream to Ignite the combustible mixture In the 
cold fan stream. The present research compares the capability for soft Igni- 
tion and the ability of different configurations of piloted flameholders to 
Ignite a cold fan stream. 


SYMBOLS 

C pressure loss coefficient 
F area 

G air-mass flow rate 

a 

G f fuel flow rate 

h gap height of V-gutter or peripheral length of test section 

M Mach number 

2 

P pressure, kg/cm 

3 

P* stagnant pressure, kg/cm 

2 

AP pressure drop or pressure oscillation value, kg/cm 
T temperature, °C 

RECEDING PACE ,£LANX NOT FILMED 
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T* stagnant temperature, °C 
V velocity, m/sec 
AFR air-fuel ratio 
Subscripts : 

0 upstream of flow orifice 

1 downstream of flow orifice or upstream of air heater 

2 Inlet of test section 

3 outlet of test section 
a air or ambient 

ave average 


APPARATUS 

The test apparatus Is shown In figure 1. The airflow from the compressor 
Is divided Into two streams: one stream goes to ejector 7 to create a low- 

pressure exhaust condition for the test section; the other stream goes through 
valve 1, flow-orifice meter 2, air heater 3, control valve 4, and plenum cham- 
ber 5 Into test section 6. After leaving the test section the combustion 
products are cooled by water Injection. The hot gas Is exhausted Into the 
ambient through the ejector. 

The test rig assembly Is shown In figure 2. It Is a rectangular duct 
300 mm by 175 mm. The piloted flameholder Is mounted Inside the test section 
by bolts. Flame extinction Is recorded through a quartz window by a flame- 
detection device. Flame extinction can also be observed visually. The Igni- 
tion In the primary zone Is provided by a high-energy Ignitor of 12 J. The 
combustible mixture In the secondary zone Is Ignited by flame propagation from 
the primary zone. The Inlet flow Is measured by the static pressure tap, the 
stagnant pressure tube, and the chromel-alumel thermocouple. The outlet tem- 
perature Is monitored by a platlnum/platlnum-rhodlum thermocouple. Outlet 
stagnant pressure Is measured by three water-cooled stagnant pressure probes; 
each has three measuring points. The transient pressure Is measured by an 
Induction-type pressure transducer (CyG-1, 0 to 1.0 kg/cm^) with an LR-1 
carrier wave amplifier and an Sc-60 ultraviolet oscilloscope. The fuel flow 
rate Is measured by a turbine-type flowmeter with a digital frequency 
Indicator. 

The fuel Is Injected through a spray bar. There Is one spray bar In the 
primary zone. The total length of the bar Is 280 mm. There are six pairs of 
holes (0.4 mm In dlam) on the bar at an Intersection angle of 60°. In the 
secondary zone there are six spray bars. The two bars In the center have four 
pairs of holes (0.35 mm In dlam) at an Intersection angle of 30°; the other 
four spray bars have four holes (0.35 mm In dlam) each. Fuel Injection In the 
secondary zone Is In the counterflow direction. The test configurations of 
piloted flameholders are shown In figure 3. 
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TEST METHOD 


The test conditions were as follows: 

T 2 60 to 70 °C or 500 °C 
M 2 0.23 

P 2 0.25 to 1.2 kg/cm 2 

The test procedures were as follows: First, the air supply was adjusted to 

maintain the Inlet conditions of P 2 , M 2 , and T 2 as specified. After the 
primary zone had been Ignited the ejector air supply was adjusted to maintain 
P 2 , M 2 , and T 2 constant. Then the fuel flow In the primary zone Gf p was 
Increased until fuel-rich flameout occurred. After that, the primary zone was 
reignited at the specified Inlet condition and P 2 , M 2 , and T 2 were kept 
constant. Finally, the fuel flow was reduced until there was fuel-lean flame 
extinction. Before ignition, the pressure drop between the Inlet and outlet 
was measured. During ignition the peak pressure, the fuel-air mixture In the 
secondary zone, and the pressure oscillation were also recorded. 


RESULTS AND DISCUSSION 
Peak Pressure Measurements 

As shown In figure 4, the peak pressure rise during Ignition is usually 
lower than 10 percent of the Inlet pressure (AP 2 /P 2 < 10 percent). The 
peak pressure rise in the augmentor is normally related to the amount of fuel 
Injected (as In these tests). Since the Ignition is realized in a small vol- 
ume, the amount of fuel necessary for Ignition Is small. Thus, the peak pres- 
sure is lower than that for normal augmentor Ignition. This proves that using 
a piloted flameholder makes It easier to obtain a soft ignition in a turbofan 
augmentor. 


Pressure Oscillation Measurements 

Pressure oscillation measurements were performed under the same conditions 
as the peak pressure measurements. The amount of fuel Injected in the primary 
zone has a very strong effect on the pressure oscillation in the augmentor 
during ignition In the secondary zone (figs. 5(a) and (b)). It has been shown 
that, with the same amount of fuel Injected in the secondary zone, the pressure 
oscillation during the onset of secondary combustion varies with the amount of 
fuel Injected In the primary zone. It Is clear that there Is an optimum fuel 
injection in the primary zone which minimizes the pressure oscillation during 
secondary zone ignition. 


Flow Resistance 

Figure 6 shows the variation of the pressure loss coefficient with the 
Reynolds number Re for different configurations. When Re Increases, the 
pressure loss coefficient C decreases slightly. The pressure loss coeffi- 
cient C also decreases with the increase of gap height h. The pressure loss 
coefficient C for a V-gutter is slightly lower than that for a semi spherical 
shape. 


193 



Effect of Inlet Flow Distortion 


The inlet flow distortion is shown in figure 7. Figure 8 shows the effect 
of inlet flow distortion on flame stabilization. When the inlet flow is dis- 
torted, the range of flame stabilization is narrowed, and the peak value of the 
stabilization parameter is also reduced. For instance, at P 2 = 0.25 kg/cm 2 
when there is no inlet distortion, it is possible to ignite the secondary zone 
in 3 sec; but, when there is inlet distortion (even at P 2 = 0.3 kg/cm 2 ), the 
range of flame stabilization is very narrow and the ignitability in the sec- 
ondary zone is very weak. (It is impossible to ignite the fuel-air mixture in 
the secondary zone.) 

Effect of Secondary Air Entry Holes on Stabilization 
of the SJ410-4A Flameholder 

When 54 secondary holes (8 mm in diam) are opened, the range of flame 
stabilization is very narrow (fig. 9). When these holes are blocked, the 
stabilization range is broader. The secondary air entry holes also reduce the 
capability of igniting the fuel-air mixture in the secondary zone. It seems 
that the secondary air dilutes the primary combustion products and reduces the 
temperature. 


Effect of Gap Height 

Figure 10(a) and (b) show the effect of gap height h on flame stabili- 
zation for the SJ410-4B and SJ410-3A configurations. The range of flame 
stabilization decreases with the Increase of the gap height h, and the Ignit- 
ability also decreases. From this test we found that, for a specified con- 
figuration, there is an optimum gap height. For Instance, when the width of 
the V gutter is 70 mm (SJ410-3A), the optimum gap height is 3 to 5 mm, the 
width of V-gutter is 40 mm (SJ410-4B), and the optimum gap height is 2 to 3 mm. 
Under these conditions, the range of ignition and flame stabilization and the 
ignition capability in the secondary zone are nearly the same for the two 
conf igurations . 


Effect of Bluff Body Shape on Flame Stabilization 

Figure 11(a) shows the flame stabilization range for the SJ410-3A con- 
figuration under ambient air temperature with a different bluff body shape. 

The range of flame stabilization of the V-gutter bluff body is wider than that 
of the semispherical bluff body. But as shown in figure 11(b), at t£ = 500 °C 
the ranges of flame stabilization for both bluff bodies are nearly the same. 
This can be explained as follows. At high inlet temperature, the. vaporization 
rate of the fuel spray (before entering the primary zone) is high for both 
f lameholders ; therefore, the gaseous fuel-air ratios of the mixtures in the 
primary zone are nearly the same. But at low temperature, the fuel vapori- 
zati on r d 1 0 on tho y - 5 u 1 1 0 r s u r f d c 0 is hiQhcr thdn thdt on th 6 sorni-sphoricdl 
body, and the fuel-air ratio in the primary zone of the V-gutter is higher than 
that of the semi-spherical body. Thus, the V-gutter has a wider range of flame 
stabi 1 i zation. 

Effect of Air Entry Holes in V Gutter Wall on Flame Stabilization 


Figure 12 shows the range of flame stabilization for the SJ410-3B con- 
figuration with different air entry holes in the V-gutter wall. For a given 
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opening area, fewer holes with larger diameters offer better flame stabili- 
zation than smaller holes. A comparison of figure 10(b) with figure 12 shows 
that the recirculation zone created by the gap provides better flame stabili- 
zation than that created by circular holes. 


CONCLUSIONS 

These tests led to the following conclusions. 

(1) The use of a piloted flameholder In the turbofan augmentor may 
minimize the peak pressure rise during Ignition. At the present experimental 
conditions, AP/p£ Is less than 10 percent; therefore, the use of a piloted 
flameholder Is a good method to realize soft Ignition. 

(2) The geometry of the piloted flameholder and the amount of fuel 
Injected Into the flameholder have a strong effect on the pressure oscillation 
during Ignition of the fuel-air mixture In the secondary zone. 

(3) Compared with the V-gutter flameholder with holes In Its wall, the 
V-gutter flameholder without holes not only has the advantages such as simple 
structure and good rigidity but offers a wide combustion stability limit and a 
high capability of Igniting the fuel-air mixture of the secondary zone. 
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Air supply 



Figure 1. - Test apparatus. Valve, 1; orifice flowmeter, 
2; air heater, 3; control valve, 4; plenum chamber, 

5; test section, 6; ejector, 7; valve, 8. 



Figure 2. - Test rig assembly. High energy Ignitor, 1; quartz 
window, 2; thermocouple, 3; pressure probe, 4; pressure 
transducer, 5; test configuration, 6; fuel spray, 7; fuel 
spray bar, 8; thermocouple, 9; pressure probe, 10; pressure 
tap, 11 . 
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1 1 7 7 7 1 7 7 , T 


(a) Configuration SJ410-4A; h = 4 mm. 



(b) Configuration SJ410-4B; h = 2, 3, 4, and 6 mm. 

Figure 3. - Test configurations of piloted f lameholders . 
(All dimensions in mm.) 




Type 


Number 

Diameter 

A 

■ 


1 2mm 

B 

1 4 

76 

8 cam 

C 

7 

304 

4 mm 


(e) Air-entry holes on the V-gutter. 


Figure 3. - Concluded. 


198 













Type 

A 

B 

C 

D 

E 

P 2 kg/cm 2 

1.2 

o.e 

o.e 

0.4 

0.3 

kg/ca 2 

0.057 

0.047 

0.051 

0.0085 

0.0085 

— 1 

^ P 2 /f 2 * 

4.8 

5.9 

2.5 

2. 1 

2.8 



Figure 4. - Pressure oscillation during ignition. 
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(a) Pressure oscillation during combustion start- 
ing in the secondary zone. 




















E 


Figure 10. - Effect of gap height on flame stabilization. 










Figure 12. - Effect of air-entry holes in 
V-gutter wall on flame stabilzation. 
Configuration SJ410-3B; T 2 * = 500 °C, 

M 2 = 0.23. 
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THE CHEMICAL SHOCK TUBE AS A TOOL FOR STUDYING HIGH- 
TEMPERATURE CHEMICAL KINETICS 

Theodore A. Brabbs 
NASA Lewis Research Center 
Cleveland, Ohio 


Although the combustion of hydrocarbons Is our primary source of energy 
today, the chemical reactions, or pathway, by which even the simplest hydro- 
carbon reacts with atmospheric oxygen to form CO 2 and water may not always 
be known. Furthermore, even when the reaction pathway Is known, the reaction 
rates are always under discussion. The shock tube has been an Important and 
unique tool for building a data base of reaction rates Important In the 
combustion of hydrocarbon fuels. 

The ability of a shock wave to bring the gas sample to reaction conditions 
rapidly and homogeneously makes shock-tube studies of reaction kinetics 
extremely attractive. In addition to the control and uniformity of reaction 
conditions achieved with shock-wave methods, shock compression can produce gas 
temperatures far In excess of those In conventional reactors. Argon can be 
heated to well over 10 000 K, and temperatures around 5000 K are easily 
obtained with conventional shock-tube techniques. Experiments have proven the 
validity of shock-wave theory; thus, reaction temperatures and pressures can 
be calculated from a measurement of the Incident shock velocity. 

This report presents a description of the chemical shock tube and auxil- 
iary equipment and of two examples of kinetic experiments conducted In a shock 
tube. Shock-wave theory and shock tubes In general are discussed In 
references 1 and 2 and thus are not Included In this report. 


EXPERIMENTAL APPARATUS AND PROCEDURES 

The chemical shock tube used was a single piece of square stainless-steel 
tubing (5.7 m long, 6.4 cm on a side, 1.3 cm thick walls). The entire length 
of the tube was ground to a constant Inside dimension and then honed to a 
highly polished finish. Stations for shock-wave detectors were located at 
15-cm Intervals In the downstream portion of the tube. Two 2.5-cm calcium 
fluoride windows were located 25 cm from the end of the tube. A quartz pres- 
sure transducer was located at a position where It would provide an accurate 
Indication of the time at which a shock wave arrived at the center of the win- 
dows. The assembled tube could be evacuated to a pressure of about 1 ym and 
had a leak rate of less than 0.2 ym/mln. A liquid nitrogen cold trap In the 
vacuum line guarded against the back-migration of pump oil Into the shock tube. 

The temperature and pressure behind the Incident and reflected shock waves 
were calculated from the measured velocity of the Incident shock wave with the 
shock program of Gordon and McBride (ref. 3). Since the reaction temperature 
had to be calculated. It was very Important to measure the Initial temperature 
and the shock velocity as accurately as possible. The distances between probe 
locations were accurately measured to four significant figures. The shock-tube 
wall temperature was measured to ±0.5 C, and It was assumed that the test mix- 
ture rapidly achieved this temperature. Thln-fllm resistance gauges (1.3 cm 
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dlam) were used to detect the passage of the Incident shock wave. The ampli- 
fied output from these gauges was displayed on a digital oscilloscope (see 
fig. 1). Data points were recorded every 200 nsec. Since with this scope the 
scale can be magnified to allow viewing of each point, we were able to measure 
the Incident shock velocity to better than ±0.2 psec (0.1 percent). 

All data were taken at an Initial test gas pressure of 8.00 kPa (60 torr) 
and a driver pressure of 792.7 kPa (115 psla). Soft aluminum diaphragms sepa- 
rated the high-pressure driver gas from the test gas. These diaphragms were 
pressurized and then pierced by a gas-activated piston. The shock strength 
was varied by adding small amounts of argon to the helium driver gas to change 
Its molecular weight. Figure 2 shows the variation of shock velocity with the 
reciprocal of the molecular weight of the driver gas. The constant pressure 
ratio across the diaphragm Insured that the opening process for the diaphragm 
was about the same for all tests. 

Gas mixtures were prepared by the method of partial pressures In 34.4- 
1 1 ter stainless-steel tanks. The gases used to prepare the mixtures had stated 
purities of 99.99-percent methane, 99.98-percent oxygen, and 99.998-percent 
argon. The partial pressure of methane and oxygen could be determined to 
better than 0.3 mm of Hg, and the total pressure of the tank to better than 
1.38 kPa (0.2 psla). All samples were prepared to a total pressure 
of 413.6 kPa (60 psla). 

Radiation at 3700 A was monitored for each experiment and was displayed 
on the oscilloscope with the pressure history. The radiation was detected by 
a photomultiplier tube whose field of view was restricted by two slits (0.75 
and 1.5 mm wide) located 10 and 23 cm from the center of the shock tube. With 
this arrangement we could view a 3-mm slice of the hot gas In the center of 
the shock tube at the same location as the pressure transducer. An Interfer- 
ence filter centered at 3696.5 A (with a half band width of 115 A) was used to 
Isolate the 3700 A radiation. The filter transmitted 51 percent of the Inci- 
dent radiation. 

Examples of kinetic experiments conducted In a shock tube are presented 
In the next two sections. 


SHOCK TUBE MEASUREMENTS OF SPECIFIC REACTION RATES IN THE BRANCHED-CHAIN 

H 2 -C0-0 2 SYSTEM 

When a mixture containing hydrogen, carbon monoxide, and oxygen Is sub- 
jected to a temperature and pressure pulse In a shock tube, small concentra- 
tions of atoms and free radicals are first formed by processes such as 

CO + 0 2 - » C0 2 + 0 

These concentrations may then grow exponentially via the well-known branched- 
chaln scheme 


OH + H, 


h 2 o 


+ H 


I 
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OH t 0 


II 



III 

IV 


OH + CO — ^—C0 2 f H V 

Theory shows (refs. 4 and 5) that the oxygen atom concentration grows as 
[0] sr A exp xt, (except very early In the reaction). Where X Is the posi- 
tive root of the cubic equation 

3 2 

X +• (v-| * v 2 4- w 3 + + Vg)X 4- [( V] + v & )(v 3 4- « 4 ) 4- u-jV^X 


- v 3 ( v i + v 5^ 2v 2 " v 4> = 0 (!) 

where v-| .. k-|[H 2 ], v 2 = ^2^2^’ v 3 ~ ^3^2^’ ~ ^ 4 ^ 2 ^^’ w 5 — kg[C0]. 

Thus, the exponential growth constant x depends on the gas composition and 
the rate constants of reactions I to V. This paper reports measurements of 
mixtures chosen to permit determination of the rates of reactions I, II, III, 
and V. Mixtures were selected by analyzing equation (1). 

Growth constants were obtained by measuring the blue carbon monoxide 
flame-band emission behind Incident shocks. The Intensity of this radiation 
Is proportional to the product of carbon monoxide and oxygen atom concentra- 
tions (ref. 6), and, since very little carbon monoxide Is consumed, the light 
monitors the Increase of oxygen atom concentration with time. A typical emis- 
sion record Is presented In figure 3. 

Exponential growth constants were obtained from plots of the logarithm of 
observed light Intensity as a function of gas time (see fig. 4); the relation 
between gas and laboratory times was obtained from the computer calculations. 


SELECTION OF GAS MIXTURES 

Let us turn now to the rationale used to select gas mixtures by analyzing 
equation (1). To begin with, under our experimental conditions V 4 Is 
generally small In comparison with other v; therefore. It can be neglected 
for purposes of a qualitative discussion. Secondly, X turns out to be a small 
positive root - of the order of the smaller . y values and very small compared 
with the larger v values. Thus, we neglect x^ In comparison with the 
other terms and rewrite equation ( 1 ) as 

2 

[(”l + V 5 ) 4- «2 + v 3^ + f ~ 2 w 2 u 3^ w 1 + v 5^ ( 2 ) 
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If the amount of hydrogen In a mixture Is large In comparison to the amount of 
oxygen, v-\ and V 3 are large and the term Involving X 2 may be neglected; 

In this event, 


X = 2v 2 (3) 

On the other hand. If only a trace of hydrogen Is present; *>3 Is small, the 
term Involving X may be neglected, and 


2 e 2 v 2 v 3 (v 1 * V 

v 2 + (v l + V 


(4) 


If we choose a mixture with a large amount of carbon monoxide, V 5 Is large and 

X ~y2u 2 v 3 (5) 

If there Is a large amount of oxygen, v 2 Is large and 

x ~y2v 3 (v 1 ~ v 5 ) 

X-^ for [H 2 ] > [CO] ( 6 ) 

X for [CO] > [H 2 ] 

This, then, outlines a strategy for obtaining rates of reactions I, II, 
III, and IV. First, a mixture rich In hydrogen Is used to determine k 2 . 

Next, with k 2 known, a mixture with a trace of hydrogen and rich In carbon 
monoxide Is used to determine k 3 . Finally, with k 3 known, mixtures with 
excess oxygen and varying proportions of hydrogen and carbon monoxide are used 
to Isolate k-] and k 5 . 


The foregoing discussion Indicates a qualitative procedure for selecting 
gas mixtures. However, we also need a quantitative measure of the sensitivity 
of the growth constant for a particular mixture to the various rate constants. 
For example, we know that to Isolate k 3 we need a mixture with a trace of 
hydrogen, a small amount of oxygen, and an excess of carbon monoxide. But, for 
a candidate composition, Is the hydrogen concentration small enough, and Is 
carbon monoxide sufficiently In excess? A quantitative measure of the sensi- 
tivity to the various rates can be obtained by logarithmic partial differentia- 
tion of equation (1). For example, 



1 

21 


X + 2v 3 X + ( v-| <■ Vg)(vg - V 4 ) + 

- 

3X + 2(v n + v 0 + v Q + Vy, + v K ) + (u-, + + v A ) + 


3 4 


(7) 


Equation (7) shows that 
yiTJ (a somewhat weaker 


the growth constant depends on something less than 
dependence than that suggested by equation (5)). 
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The mixture compositions selected are shown in table I along with the 
sensitivities calculated for a pressure of 1 atmosphere and the temperatures 
near the midpoints of the ranges of the experimental data. These mixtures 
were chosen on the basis of ( 1 ) high sensitivity to the rate constant being 
determined, ( 2 ) minimized sensitivity to the other rate constants, ( 3 ) growth 
constants in a convenient range for measurement ( 5 x 1 0 3 to 7 xT 0 4 sec -1 ), and 
(4) sufficient content of minor, rate- limiting constituents to permit accurate 
mixture preparation. 

It was necessary to Include in each mixture sufficient infrared-active gas 
so that the arrival of the contact surface could be detected by the cessation 
of infrared emission. (The position of the contact surface between the driver 
and the driven gas was needed for the boundary layer analysis.) Mixtures 2, 3, 
and 5 have substantial carbon monoxide concentrations; carbon dioxide was added 
to mixture 1 to achieve adequate Infrared activity. Carbon dioxide was also 
added to mixtures 3 and 5 to ensure vibrational relaxation of the carbon mon- 
oxide. In preliminary experiments on a composition similar to mixture 3 (both 
without C0 2 ) carbon monoxide was not relaxed and rate constants gave an acti- 
vation energy of 15 kcal for reaction 3; this is much higher than the values 
shown in reference 7. Results from mixture 3, with added C0 2 , show that the 
preliminary data were too low (particularly at the lower temperatures). 


CALCULATION OF RATE CONSTANTS FROM EXPONENTIAL GROWTH CONSTANTS 

Rate constants were obtained from the experimental growth constants by 
means of equation ( 1 ), which was rearranged and solved for the rate constant 
being sought from the mixture in question: k 2 from mixture 2 , k 3 from mixture 
3, etc. Mixture 2 was studied first, because the growth constants for this 
mixture depend almost exclusively on k 2 and are only slightly affected by 
the rates of the other reactions (table 1 ). Values of «4 were calculated 
from the k 4 suggested by the Leeds Group (ref. 8 ). Trial values of k-| , k 3 , 
and k 5 were also taken from the Leeds recommendations (refs. 7 to 9), but 
final values were taken from our own determinations. This Involves an itera- 
tive procedure since our values of k 3 depend, in turn, on our values of k 2 
and, to a lesser extent, of k 5 . Two iterations sufficed to establish k 2 . 

Next, k 3 was obtained by analyzing the growth constants for mixture 3, 
then k] was obtained from mixture 1 , and finally k 5 was obtained from 
mixture 5. Three Iterations around the k 3 k 3 - k 5 loop were required. 

A least-squares fit to the Arrhenius equation was made for each set of 
rate constants; these equations were then used in subsequent calculations. 

The experimental results are presented in the order in which the data 
were obtained - first k 2 , then k 3 , k] , and k 5 . This is also the order of 
decreasing precision. 

Rate constants for reactions 11 and III are plotted as functions of recip- 
rocal temperature in figure 5. The least squares lines through the data are 

k 2 = 1.25xlCP 4 exp(-16.3 kcal/RT) cm^ mole - ! sec - ! 

k 3 = 2.96x1()!3 exp(-9.8 kcal/RT) cm^ mole - ! sec - ! 
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Figure 6 contains the data for reactions I and V plotted as functions of the 
reciprocal temperature. The least-squares lines through the data are 

k-j = 1.9xl0 13 exp(-4.8 kcal/RT) cm 3 mole -1 sec'] 

k 5 = l.OxlO 12 exp( -3 .7 kcal/Rl ) cm 3 mole' 1 sec' 1 

In this paper we showed how to make systemmatlc measurements so as to 
Isolate the rates of Individual blmolecular reaction steps In the complex of 
six reactions which describe the Ignition behavior of the hydrogen-carbon 
monoxide-oxygen system. 

We have obtained what we believe are quite accurate measurements of the 
rates of the reactions H + 0 2 — ►OH + 0 and 0 + H 2 -»-0H + H. Proper accounting 
of the effects of boundary layer growth on the pressure, temperature, and 
residence time behind the shock wave has been essential. Our measurements of 
the rates of OH + H 2 — ►H 2 0 + H and OH + CO -**C0 2 + H were less precise but 
were perhaps the most direct determinations of these rates for temperatures 
above 1100 K (previous high-temperature determinations of these rates have been 
based on flame sampling). Finally, we have estimated the rate of the Initia- 
tion reaction CO + 0 2 — ►C0 2 + 0. 


METHANE OXIDATION BEHIND REFLECTED SHOCK WAVES, IGNITION DELAY TIMES 
MEASURED BY PRESSURE AND FLAME-BAND EMISSION 

Large computers have enabled us to assemble kinetic mechanisms for model- 
ing the combustion of many fuels. These mechanisms may have as many as 150 
reactions (variables) and In some cases the only experimental parameter matched 
has been the Ignition delay times obtained from shock tubes. Many studies 
have been conducted on the Ignition of methane-oxygen gas mixtures In shock 
tubes. However, In some cases the data are In error because of probe location, 
and In other cases the data are badly scattered because of poor measurements 
of the incident shock velocity or a very noisy pressure trace. 

The delay time measured for a shock-heated hydrocarbon-oxygen mixture Is 
known to be a function of Initiation reactions, of reactions of radicals and 
molecular oxygen with the fuel, and of branching reactions. These branching 
reactions Increase the radical concentrations to a level where Ignition will 
occur. During Ignition, there Is rapid depletion of the primary fuel, very 
high radical concentrations, and an exponential rise In temperature and pres- 
sure. Usually delay times are measured somewhere In this region by techniques 
which follow some physical process such as change In pressure or In the appear- 
ance of some emitting species. Although the time difference between the start 
and completion of the pressure rise Is small, one would like to determine the 
time of the appearance of the signal as accurately as possible. This time, 
which Is the start of the Ignition process, allows one to separate the fast 
Ignition mechanism from the kinetics of the primary fuel. 

The next section presents Ignition delay data for the oxidation of 
methane. Data were measured with a quiet pressure transducer for the tempera- 
ture range of 1500 to 1920 K. The quiet pressure trace allowed us to determine 
the Initial rise In the Ignition pressure. The delay times determined from 
this Initial pressure rise were compared with the time of appearance of the 


212 


radiation at 3700 A. In addition, the methane and carbon dioxide concentra- 
tions at the time of Ignition will be measured. These parameters, which must 
be duplicated by any kinetic model, will be used to test the correctness of 
published kinetic models for the combustion of methane. 


RESULTS AND DISCUSSION 
Ignition Delay Times and Probe Location 

Frenklach (ref. 10) measured Ignition delay times at two locations In a 
shock tube and concluded that the location of the probe must be taken Into 
consideration when selecting delay-time data to be modeled by a kinetic mecha- 
nism. Since optical access was necessary to monitor the Infrared and ultra- 
violet radiation behind the reflected shock wave. It was Important to know the 
error Involved In measuring the delay time 7 mm from the reflecting surface. 
Therefore, two pressure transducers were located 7 and 83 mm from the reflect- 
ing surface. These allowed us to simultaneously measure the delay time for 
two locations In the same experiment (see fig. 7). An ethane-oxygen-argon 
mixture (1.27, 5.05, and 93.68 percent, respectively) was used for this series 
of tests. For these experiments, delay time was defined as the time at which 
the Ignition pressure exceeded P 5 by 10 percent. This eliminated the 
guess-work usually Involved In determining the point at which the pressure 
departed from the flat portion of the trace. Figure 8 Is a plot of the delay 
times for the two locations as a function of 1/T. A curve drawn through the 
data taken at the 83-mm position appears to show a constant difference In the 
time between the two locations. 

A clearer understanding of this phenomenon can be obtained from figure 9, 
which Is a time-distance (t-x) plot of the trajectory of the Incident and 
reflected shock waves near the end of the shock tube. As can be seen, the 
time at which shock heating occurs at a given location along the tube Is a 
function of the reflected shock velocity and the distance from the end wall. 
Thus, Ignition would be expected to occur at the end wall first, since the gas 
there has been heated much longer. When Ignition occurs, the pressure rise 
causes a disturbance which propagates down the tube. The mlnlnum propagation 
velocity of this disturbance Is the speed of sound In the shocked gas. Vermeer 
(ref. 11) took schlleren photographs of Ignition behind the reflected shock 
wave, and they clearly showed that the disturbance propagated at a velocity 
much greater than the reflected shock wave. 

The data In figure 8 were used to calculate the propagation velocity 
between the two probes. This velocity had an average value of 1660 m/sec, 
which Is a Mach number of about 2.5. This velocity Is about four times that 
of the reflected shock wave. Using this velocity one can estimate that the 
delay times measured at the 7-mm position would be too short by about 10 ysec. 
This error would not be Important unless the measured delay times were near 
100 ysec. No corrections have been made In the present data, but times less 
than 100 ysec are not considered reliable. 

Another disturbing feature of the pressure traces Is the amplitude of the 
Ignition pressure. Unless something changes, one expects the Ignition pressure 
to be the same at all locations In the tube. However, the Ignition pressure 
for probe 2 Is much larger than that at probe 1 (fig. 4). The ratio of the 
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Ignition pressure to the reflected shock pressure P5 for the two probe 
positions Is plotted as a function of the delay time at probe 1 In figure 10. 
Note that this ratio Is nearly constant at probe 1 but Increases with decreas- 
ing delay times for probe 2. This behavior suggests that the gas properties 
at Ignition are not the same at the two locations. 


Ignition Delay Times 


Ignition delay times were measured behind reflected shock waves for the 
mixtures shown In table II. The pressure history and the carbon monoxide 
flame-band emission at 3700 A were used to determine the onset of Ignition. 

The delay times measured by both pressure and flame-band emission are recorded 
In table III. These delay times t may be correlated by an empirical 
relation of the form 


t = A exp(B/RT)[CH 4 ] a [0 2 ] b [Ar] c (8) 

LI f shl tz (ref. 12) and Tsubol (ref. 13) both found an overall pressure depend- 
ence of -0.7 for a + b * c. Plots of tP 0 * 7 as a function of 1/T for each 
mixture are shown In figure 8. A least-squares fit to the data yielded the 
following equations: 


= 0.5 

,p 0 - 7 

= 9 . 59x1 0 -5 

exp( 51 .93 x10 3 /RT) 

. 0.7 
ysec-atm 

(9) 

= 1.0 

tP 0.7 

= 2.85xl0~ 4 5 

exp(58.83x!0 3 /RT) 

. 0.7 
ysec-atm 

(10) 

= 2.0 

,P 0 - 7 

= 3.65X10" 5 

exp(58.74x!0 3 /RT) 

. 0.7 
ysec-atm 

(11) 


A correlation of all of the data In the form used by LI f shl tz , 
t [ CH 4 ]-° * 33 [ O2 ] 7 - 03 , Is shown In figure 11. The activation energies for the 
present data were larger than those found by Llfshltz and appear to have 
depended upon the concentration of oxygen. The mixtures of <p =1.0 and 2.0 
had 4-percent oxygen and the same activation energy, while the mixture for 
9 = 0.5 had 8-percent oxygen and a measurably lower activation energy. This 
is the same behavior reported by Tsubol. 
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TABLE 1. - MIXTURE COMPOSITIONS AND 
GROWTH CONSTANT SENSITIVITIES 


Constituents 

Mixture 


1 

2 

3 

5 


Reaction 


OH + H- 

H + 0- 

0 * H- 

OH CO- 


H 2 0 + H 

OH + 0 

OH + H 

C0 2 + H 


Composition, percent 

h 2 

0 .21 

5 

0.1046 

0.1035 

CO 

.11 

6 

10.0 

6.01 

°2 

10.0 

0.5 

.503 

10.0 

C0 2 

5.0 

- - 

4.99 

5.0 

Sensitivities 

Mixture 


1 

2 

3 

5 


Reaction 


OH + H£ 

H + 0- 

0 + H- 

OH CO- 


H 2 0 + H 

OH + 0 

OH + H 

C0 2 + H 

(8 In x) /( 8 Tn v*| ) 

0.34 

0.01 

0.00 

0.07 

(8 Tn X) /( 8 Tn v^) 

.33 

1 .00 

.64 

.21 

(8 In X) /( a In v 3 ) 

.48 

.06 

.39 

.49 

(8 In X)/( 8 In v^) 

-.17 

-.07 

-.06 

-.06 

(8 In X) /( 8 In v c ) 
□ 

.02 

.00 

.04 

.29 
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TABLE II. - COMPOSITION OF GAS 
MIXTURES AND TEMPERATURE 
RISE DUE TO IGNITION 


Equivalence 

ratio, 

«p 

ch 4 

o 2 

Ar 

Temperature 
rise due to 
Ignition, 

K 

Composition, 

percent 

0.50 

19 

8 

90 

894 

1 

■9 

n 

94 

882 

2 

H 

a 

92 

729 


TABLE III. - IGNITION DELAY TIMES MEASURED BY 

PRESSURE AND FLAME BAND EMISSION 

(a) Equivalence ratio, <p, 0.5; methane, 2 percent; 
oxygen, 8 percent; argon, 90 percent 


Temperature, 

K 

VT, 

K" 1 

Pressure, 

atm 

Delay times, tisec 

Pressure 

Emission 

1778 

5. 624x1 0~ 4 

3.464 

99 

110 

1763 

5.672 

3.418 

102 

110 

1752 

5.708 

3.371 

112 

125 

1730 

5.780 

3.318 

159 

167 

1695 

5.900 

3.199 

201 

208 

1687 

5.928 

3.191 

208 

209 

1657 

6.035 

3.072 

323 

322 

1620 

6.173 

2.979 

508 

512 

1610 

6.211 

2.950 

533 

542 

1605 

6.231 

2.918 

594 

597 

1601 

6.246 

2.922 

672 

672 

1579 

6.333 

2.857 

663 

658 

1576 

6.345 

2.835 

742 

738 

1555 

6.431 

2.773 

850 

858 

1554 

6.435 

2.786 

958 

964 

1531 

6.532 

2.704 

1144 ! 

1145 

1499 

6.671 

2.613 

1655 

1649 
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TABLE III. - Concluded. 


(b) Equivalence ratio, tp, 1.0; methane, 2 percent; 
oxygen, 4 percent; argon, 90 percent 


Temperature, 

K 

1/1, 

K-l 

Pressure, 

atm 

Delay times, ysec 

Pressure 

Emission 

1915 

5 .222x1 O' ^ 

3.680 

57 

59 

1866 

5.359 

3.536 

84 

78 

1802 

5.549 

3.354 

171 

166 

1780 

5.618 

3.279 

201 

219 

1775 

5.638 

3.278 

185 

192 

1746 

5.727 

3.183 

295 

295 

1732 

5.774 

3.147 

383 

383 

1725 

5.797 

3.122 

397 

396 

1710 

5.848 

3.082 

427 

414 

1707 

5.858 

3.077 

531 

559 

1689 

5.921 

3.023 

525 

522 

1689 

5.921 

3.032 

589 

589 

1655 

6.042 

2.932 

794 

775 

1645 

6.079 

2.900 

867 

861 

1635 

6.116 

2.876 

967 

973 

1629 

6.139 

2.855 

1072 

1075 

1615 

6.192 

2.822 . 

1163 

1144 

1593 

6.278 

2.761 

1448 

1471 


(c) Mixture ratio, <p, 2.0; methane, 4 percent; 
oxygen, 4 percent; argon, 92 percent 


Temperature, 

K 

1/1 , 
K-l 

Pressure, 

atm 

Delay times, ysec 

Pressure 

Emission 

1922 

5 . 203x1 0" ^ 


67 

62 

1856 

5.388 

3.797 

116 

110 

1789 

5.590 

3.583 

198 

201 

1770 

5.650 

3.524 

249 

271 

1763 

5.672 

3.485 

318 

304 

1730 

5.780 

3.381 

369 

368 

1715 

5.831 

3.351 

509 

501 

1712 

5.841 

3.326 

489 

486 

1691 

5.914 

3 . 2bU 

690 

682 

1681 

5.949 

3.228 

765 

767 

1677 

5.963 

3.217 

795 

797 

1650 

6.061 

3.134 

1102 

1096 

1627 

6.146 

3.063 

1214 

1221 

1621 

6.169 

3.043 

1319 

1321 

1602 

6.242 

2.984 

1534 

1557 
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ORIGINAL PAGE 

OP POOR QUALITY 



Figure 1. - Shock velocity measurement. Output from thin-film probe 
amplifier fed to digital oscilloscope. Time per point, 200 n sec. Lower 
trace shows expansion of scale, which is helpful in determining exact 
time shock wave passed probe position. 
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ORIGINAL PAGE 5S 

OF POOR QUALITY 



Figure 7. - Delay times measured at two locations in shock tube. Probes 1 and 2 are 7 and 83 mm from 
reflecting surface, respectively. P2 is incident shock pressure and P5 is reflected shock pressure. 
Measured delay times were 712 and 574psec. 
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Figure 8. - Ignition delay times for a 1. 27 -per- 
cent CpH^, 5.05-percent Op. 93.68-percent 
Ar ( (p ■ 0. 88) mixture as recorded at two 
different locations along shock tube. 



Figure 9. - Time-distance diagram of incident and 
reflected shock trajectory and delay times t 
measured at two locations. Propagation of ignition 
wave shown by dashed line. 
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DISTANCE FROM 
END WALL, 



DELAY TIME AT PROBE 1, Msec 


Figure 10. - Ratio of ignition pressure to reflected 
shock pressure, demonstrating enhanced igni- 
tion at second probe position, particularly for 
shorter delay times. 



Figure 11. - xp°* 7 as function of 1/T for all three 
CH 4 mixtures, <p* a 5, - L0, and <0- £0. 
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